最小二乘法拟合平面（线性回归法、梯度下降、PCA法）-EW帮帮网

参考笔记：

Open3D 最小二乘拟合平面（直接求解法）【2025最新版】_python open3d已知平面方程绘制平面-CSDN博客

1.前言

2.线性回归法

4.1 最小二乘法 vs PCA 拟合平面的本质区别

4.3.3 Open3D 实现PCA法拟合

5.SVD分解法

6.知识点补充：已知平面方程，如何求法向量?

1.前言

在阅读本文之前的前置内容是：最小二乘法拟合直线，可以看我写的另一篇博客：

最小二乘法拟合直线，用线性回归法、梯度下降法实现-CSDN博客文章浏览阅读222次，点赞7次，收藏4次。对于每个数据点，预测值为最小化所有数据点的误差平方和：最终目标就是找到使总误差S最小的a、b，这可以通过S分别对a、b求偏导，并令偏导数 = 0来实现梯度下降（Gradient Descent）是一种迭代优化算法，通过不断沿损失函数负梯度方向更新参数，逐步逼近最优解。https://blog.csdn.net/m0_55908255/article/details/148018256?spm=1011.2415.3001.5331在学习了最小二乘法拟合直线之后，我们可以将最小二乘法推广至三维空间，在三维空间中为 n 组三维数据点 $\color{red}(x_1,y_1,z_1),(x_2,y_2,z_2)....(x_n,y_n,z_n)$ 拟合出一个平面 $z = {\color{red}a}x+{\color{red}b}y+{\color{red}c}$ ，使得所有数据点到平面的 垂直距离平方和 最小。如下图所示：

2.线性回归法

线性回归法也称为直接法，计算简单，可以直接推导出拟合平面方程 $z = {\color{red}a}x+{\color{red}b}y+{\color{red}c}$ 的 a、b、c，应用比较广泛

2.1 模型假设

假设我们有 n 组三维数据点 $\color{red}(x_1,y_1,z_1),(x_2,y_2,z_2)....(x_n,y_n,z_n)$ ，我们的目标是找到一个平面 $z = {\color{red}a}x+{\color{red}b}y+{\color{red}c}$ ，使得所有数据点到平面的 垂直距离平方和 最小

2.2 定义误差函数

对于每个数据点 $\color{red}(x_i,y_i,z_i)$ ，预测值为 $\hat z_i = ax_i+by_i+c$ ，则误差为：

$e_i=z_i-\hat z_i = z_i-(ax_i+by_i+c)$

目标：最小化所有数据点的误差平方和：

最终目标就是找到使总误差 S 最小的 a、b、c，这可以通过 S 分别对 a、b、c 求偏导，并令偏导数 = 0 来实现

2.3 求偏导并解方程

为了找到使总误差 S 最小的 a、b、c，可以对 a、b、c 分别求偏导，令导数为 0 ，得到方程组如下：

将方程组展开并整理为矩阵形式：

注： $\sum =\sum_{i=1}^n$

2.4 解方程

① 计算中间项：

$\sum x_i,\sum y_i,\sum y_i$
$\sum x_i^2,\sum y_i^2,\sum x_iy_i$
$\sum x_iz_i,\sum y_iz_i$

全是一些常数，直接进行计算即可

② 构建系数矩阵和常数向量：

③ 解方程组：

$\begin{bmatrix} a \\ b \\ c \end{bmatrix}=A^{-1}B$

④ 拟合平面：

$z = {\color{red}a}x+{\color{red}b}y+{\color{red}c}$

补充：解方程时也可以使用克拉默法则，这样就不需要求 A 的逆矩阵，如下：

2.5 案例演示

2.5.1 手工计算实例

假设有以下数据点：

$(1, 2, 3), (2, 1, 4), (3, 3, 7), (4, 2, 6), (5, 2, 7), (4, 6, 11)$

计算中间项：

构建方程并代入数值：

解方程组可得：

$a=0.889,\;\;b=1.162,\;\;c=0.42$

最终拟合的平面方程：

$z = {\color{red}a}x+{\color{red}b}y+{\color{red}c}=0.0889\;x +1.162\;y+0.42$

🆗，以上就是整个手工计算过程，还是比较简单的，但需要一些《线性代数》的相关基础

2.5.2 使用Python实现

我们将 2.5.1 中的手工计算例子使用 Python 来实现，并作可视化，代码如下：

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["font.sans-serif"] = ["SimHei"]  #设置字体，可以显示中文
plt.rcParams["axes.unicode_minus"] = False  # 正常显示负号

# 构造三维数据点（6个）
points = np.array([
    [1, 2, 3],
    [2, 1, 4],
    [3, 3, 7],
    [4, 2, 6],
    [5, 2, 7],
    [4, 6, 11],
])

x = points[:, 0]
y = points[:, 1]
z = points[:, 2]
n = points.shape[0]#数据点个数

# ------------------------构建系数矩阵和常数向量-----------------------------
A = np.array([[sum(x ** 2), sum(x * y), sum(x)],
              [sum(x * y), sum(y ** 2), sum(y)],
              [sum(x), sum(y), n]])

B = np.array([[sum(x * z)],
              [sum(y * z)],
              [sum(z)]])

# 方程求解,np.linalg.solve(A,B)表示求 AX = B 的解X
X = np.linalg.solve(A, B)
print('平面拟合结果为：z = %.3f * x + %.3f * y + %.3f' % (X[0], X[1], X[2]))

a,b,c = X[0],X[1],X[2]

# 可视化
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

# 设置坐标轴和标题
ax.set_xlabel('X', fontsize=12)
ax.set_ylabel('Y', fontsize=12)
ax.set_zlabel('Z', fontsize=12)
ax.set_title('最小二乘法线性回归拟合平面', fontsize=14)

min_pt = np.amin(points, axis=0)  # 获取坐标最小值
max_pt = np.amax(points, axis=0)  # 获取坐标最大值

# 绘制数据点
ax.scatter(x, y, z,#点坐标
           c='red',#点颜色
           s=60, #点大小
           label='数据点')

# 生成网格点绘制平面
xx = np.linspace(min_pt[0], max_pt[0], 100)
yy = np.linspace(min_pt[1], max_pt[1], 100)
xx, yy = np.meshgrid(xx,yy)#生成网格
zz = a * xx + b * yy + c

# 绘制平面
ax.plot_surface(
        xx, yy, zz,
        color= "blue",#平面颜色
        alpha=0.5,#透明度 0~1
        label=f'拟合平面')


plt.show()

运行结果：

3.梯度下降法

梯度下降法的实现可以参考我的另一篇博客，最小二乘法拟合直线里的讲解：

4.PCA法拟合平面

PCA 全称 Principal Components Analysis，即主成分分析技术。大家应该都听过这个名字，但对于其原理可能不太熟悉，说实话我也不懂它的原理，所以这一块我主要讲 PCA 法拟合平面的具体流程，不涉及算法原理，后面会举手工计算例子和 Python 实现的例子

4.1 最小二乘法 vs PCA 拟合平面的本质区别

在开始 PCA 拟合平面的具体过程前，必须先明确两种方法的根本差异：

	最小二乘法 (线性回归)	PCA (主成分分析)
优化目标	最小化因变量 (Z) 的预测误差	最小化数据点到平面的垂直距离
几何意义	垂直Z轴方向的误差平方和	真正的几何最短距离
算法优势	更适合预测场景：给定 xy 预测 z	PCA更能反映数据的真实空间分布

4.2 PCA算法流程

① 数据准备

给定 n 个三维数据点 $\color{red}(x_1,y_1,z_1),(x_2,y_2,z_2)....(x_n,y_n,z_n)$ ，将其拼成一个矩阵 $\color{red}X$ 表示，如下：

$X = \begin{bmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ \vdots & \vdots & \vdots \\ x_n & y_n & z_n \end{bmatrix}$

② 中心化数据

求出 $\color{red}x$ 维的均值 $\color{red}\bar{x}$ ， $\color{red}y$ 维的均值 $\color{red}\bar{y}$ ， $\color{red}z$ 维的均值 $\color{red}\bar{z}$ ，每一维的数据都对应减去该维的均值，此过程叫做中心化数据，也叫特征中心化，得到矩阵 $\color{red}X_{Centered}$ ，如下：

$X_{centered} = \begin{bmatrix} x_1-\bar{x} & y_1-\bar{y} & z_1-\bar{z} \\ x_2-\bar{x} & y_2-\bar{y} & z_2-\bar{z} \\ \vdots & \vdots & \vdots \\ x_n-\bar{x} & y_n-\bar{y} & z_n-\bar{z} \\ \end{bmatrix}$

③ 计算协方差矩阵 $\color{red}C$ ：

$\mathbf{C} = \frac{1}{n-1} \mathbf{X}_{\text{centered}}^T \mathbf{X}_{\text{centered}}=\begin{bmatrix} c_{11} & c_{12} & c_{13} \\ c_{21} & c_{22} & c_{23} \\ c_{31} & c_{32} & c_{33} \end{bmatrix}$

注： $C$ 为 $3\times 3$ 矩阵

④ 计算特征值和特征向量：

解特征方程 $C-\lambda E =0$ ，将的求出特征值按大到小排序，分别是 $\color{red}\lambda_1,\lambda_2,\lambda_3$ ，其对应的特征向量 $\color{red}\xi _1,\;\xi _2,\;\xi _3$

⑤ 选择平面法向量，得出拟合平面：

平面法向量：最小特征值 $\color{red}\lambda_3$ 对应的特征向量 $\color{red}\xi _3$ 即为拟合平面的法向量， $\color{red}\xi _3$ 的分量即为 a、b、c

拟合平面方程：

$a\;(x-\bar{x})+b\;(y-\bar{y})+c\;(z-\bar{z})=0$

其中， $\color{red}\bar{x},\;\bar{y},\;\bar{z}$ 即为前面所求的均值

🆗，以上就是整个PCA的算法流程，其实还是比较清晰的，但要搞懂原理还是有点难的，原理部分以后再说吧

4.3 案例演示

4.3.1 手工计算实例

① 数据准备

假设有以下数据点：

$(0,0,0), (1, 0,1), (0, 1, 1), (1, 1, 2)$

将其拼接成一个矩阵 $\color{red}X$ 表示，如下：

$X = \begin{bmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ x_4 & y_4 & z_4 \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & 1 & 2 \end{bmatrix}$

② 中心化数据

计算均值：

${\color{red}\bar{x}} = \frac{0 + 1 + 0 + 1}{4} = 0.5, \quad {\color{red}\bar{y}} = \frac{0 + 0 + 1 + 1}{4} = 0.5, \quad {\color{red}\bar{z}}= \frac{0 + 1 + 1 + 2}{4} = 1.0$

中心化数据， $X$ --> $\color{red}X_{centered}$ ：

${\color{red}X_{centered}} = \begin{bmatrix} 0-0.5 & 0-0.5 & 0-1 \\ 1-0.5 & 0-0.5 & 1-1 \\ 0-0.5 & 1-0.5 & 1-1 \\ 1-0.5 & 1-0.5 & 2-1 \end{bmatrix} = \begin{bmatrix} -0.5 & -0.5 & -1 \\ \;0.5 & -0.5 & 0 \\ -0.5 & \;0.5 & 0 \\ -0.5 & \;0.5 & 1 \end{bmatrix}$

③ 计算协方差矩阵 $\color{red}C$ ：

$\mathbf{C} = \frac{1}{n-1} \mathbf{X}_{\text{centered}}^T \mathbf{X}_{\text{centered}}\\ =\frac{1}{3} \begin{bmatrix} -0.5 & 0.5 & -0.5 & -0.5\\ \;-0.5 & -0.5 & 0.5 & 0.5\\ -1 & \;0 & 0 &1\\ \end{bmatrix} \begin{bmatrix} -0.5 & -0.5 & -1 \\ \;0.5 & -0.5 & 0 \\ -0.5 & \;0.5 & 0 \\ -0.5 & \;0.5 & 1 \end{bmatrix} \\ = \begin{bmatrix} \frac{1}{3} & 0 & \frac{1}{3} \\ 0 & \frac{1}{3} & \frac{1}{3} \\ \frac{1}{3} & \frac{1}{3} & \frac{2}{3} \\ \end{bmatrix}$

④ 计算特征值和特征向量：

解特征方程 $C-\lambda E =0$ ，特征值：

$\lambda_1=1,\;\lambda_2=\frac{1}{3},\;\lambda_3=0$

最小特征值 $\color{red}\lambda_3=0$ 对应的特征向量（作为拟合平面的法向量）为：

$\xi_3 = \begin{bmatrix} -1\\ -1\\ 1 \end{bmatrix}$

⑤ 选择平面法向量，得出拟合平面：

平面法向量：最小特征值 $\color{red}\lambda_3$ 对应的特征向量 $\color{red}\xi _3$ 即为拟合平面的法向量，因此：

$\color{red}a=-1,\;b=-1,\;c=1$

拟合平面方程：

所以PCA法最终拟合的平面方程为 $\color{red}z = x + y$ ，该平面的法向量为：

$\xi_3 = \begin{bmatrix} -1\\ -1\\ 1 \end{bmatrix}$

🆗，以上就是整个手工计算过程，不算难，需要一些《线性代数》的相关基础

4.3.2 使用 Python 实现

我们将 4.2.1 中的手工计算例子使用 Python 来实现，并作可视化，代码如下：

import numpy as np
import matplotlib.pyplot as plt


plt.rcParams["font.sans-serif"] = ["SimHei"]  #设置字体，可以显示中文
plt.rcParams["axes.unicode_minus"] = False  # 正常显示负号

# 1，数据准备：构造三维数据点
points = np.array([
    [0, 0, 0],
    [1, 0, 1],
    [0, 1, 1],
    [1, 1, 2],
])

# 2.中心化数据
centroid = np.mean(points, axis=0) #计算均值 [μ_x, μ_y, μ_z]
centered_points = points - centroid  #中心化数据

# 3.计算协方差矩阵C
cov_matrix = np.cov(centered_points.T)  #计算协方差矩阵C (3x3)

# 4.解特征方程，计算特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix) #解特征方程
#eigenvalues：特征向量
#eigenvactors：特征向量

# 5.选择平面向量，得出拟合平面方程
min_eigen_idx = np.argmin(eigenvalues)  # 找到最小特征值索引
normal_vector = eigenvectors[:, min_eigen_idx]  #找到最小特征值对应的特征向量，作为拟合平面的法向量[a, b, c]

a, b, c = normal_vector
d = -np.dot(normal_vector, centroid)
'''
d = - [a,b,c] dot [μ_x, μ_y, μ_z] 
  = -(a * μ_x + b * μ_y + c * μ_z)
'''
print(f"平面方程: {a:.4f}x + {b:.4f}y + {c:.4f}z + {d:.4f} = 0")

# 6.可视化
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

# 设置坐标轴和标题
ax.set_xlabel('X', fontsize=12)
ax.set_ylabel('Y', fontsize=12)
ax.set_zlabel('Z', fontsize=12)
ax.set_title('PCA法拟合平面', fontsize=14)

min_pt = np.amin(points, axis=0)  # 获取坐标最小值
max_pt = np.amax(points, axis=0)  # 获取坐标最大值

# 绘制数据点
ax.scatter(points[:,0], points[:,1], points[:,2],#点坐标(x,y,z)
           c='red',#点颜色
           s=60, #点大小
           label='数据点')

# 生成网格点绘制平面
xx = np.linspace(min_pt[0], max_pt[0], 100)
yy = np.linspace(min_pt[1], max_pt[1], 100)
xx, yy = np.meshgrid(xx,yy)#生成网格
zz = (-a*xx - b*yy - d) / c #解平面方程z

# 绘制平面
ax.plot_wireframe(
        xx, yy, zz,
        color= "blue",#平面颜色
        alpha=0.5,#透明度 0~1
        label=f'PCA拟合平面')

# 绘制法向量（从均值点出发）
normal_vector = eigenvectors[:, min_eigen_idx]  # 法向量 [a, b, c]
normalized_normal = normal_vector / np.linalg.norm(normal_vector)  # 单位化
ax.quiver(
    centroid[0], centroid[1], centroid[2],  # 起点（平面中心）
    normalized_normal[0], normalized_normal[1], normalized_normal[2],  # 方向
    color='black',#颜色
    length=0.5,  # 箭头长度（根据数据范围调整）
    arrow_length_ratio=0.1,  # 箭头头部比例
    label='法向量'
)

plt.legend()
plt.show()

运行结果：

优雅~，太优雅了😈

4.3.3 Open3D 计算点云法向量并显示

链接如下：

Open3D 法线估计(1)——计算点云法向量并显示-CSDN博客调用Open3D的radius半径限制，仅考虑距离当前点半径0.01米（或单位）内的邻近点max_nn最多邻近点个数，最多取30个最近的邻域点，避免高密度区域的冗余计算简单来说：要计算出某个点的法向量时，需要先去搜索该点的邻近点，然后根据当前和这些邻近点拟合出一个平面方程，最终该平面方程的法向量即为该点的法向量。https://blog.csdn.net/m0_55908255/article/details/148060000?sharetype=blogdetail&sharerId=148060000&sharerefer=PC&sharesource=m0_55908255&spm=1011.2480.3001.8118

5.SVD分解法

😈待更新....

6.知识点补充：已知平面方程，如何求法向量?

已知平面方程为 $\color{red}z = x + y$ ，求该平面的法向量步骤如下：

① 将方程转换为标准形式

平面的标准方程为： $\color{red}Ax+By+Cz+D=0$ ，所以将 $\color{red}z = x+y$ 转化为：

$x+y-z=0$

对比标准方程可得， $A=1,B=1,C=-1,D=0$

② 直接提取法向量

对于平面方程 $\color{red}Ax+By+Cz+D=0$ ，其法向量为：

$n=\begin{bmatrix} A \\ B \\ C \\ \end{bmatrix}$

因此平面 $\color{red}x+y-z=0$ 的法向量为：

$n=\begin{bmatrix} 1 \\ 1 \\ -1 \\ \end{bmatrix}$

③ 验证法向量的垂直性质

法向量的垂直性质：平面上任意两点的向量与法向量的点积为 0 。例如：

点 $p_1(0,0,0),p_2(1,0,1)$ 在平面上，向量 $v = p_2-p_1=(1,0,1)$
计算点积： $n\cdot v = 1*1+1*0+1*(-1)=0$ ，垂直性验证成功

④ 单位法向量（可选）

若需要单位法向量（长度为 1 ），过程如下：

计算法向量的模长：

${\color{red}\|\mathbf{n}\|} = \sqrt{1^2 + 1^2 + (-1)^2} = \sqrt{3}$

计算单位法向量：

${\color{red}\mathbf{n}_{\text{unit}}} = \frac{\mathbf{n}}{\|\mathbf{n}\|} = \begin{bmatrix} \frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{3}} \\ -\frac{1}{\sqrt{3}} \end{bmatrix}$

⑤ 最终答案

平面 $\color{red}z = x + y$ 的法向量为：

$n=\begin{bmatrix} 1 \\ 1 \\ -1 \\ \end{bmatrix}$

或单位法向量：

${\color{red}\mathbf{n}_{\text{unit}}} = \begin{bmatrix} \frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{3}} \\ -\frac{1}{\sqrt{3}} \end{bmatrix}$

最小二乘法拟合平面（线性回归法、梯度下降、PCA法）

1.前言

2.线性回归法

2.1 模型假设

2.2 定义误差函数

2.3 求偏导并解方程

2.4 解方程

2.5 案例演示

2.5.1 手工计算实例

2.5.2 使用Python实现

3.梯度下降法

4.PCA法拟合平面

4.1 最小二乘法 vs PCA 拟合平面的本质区别

4.2 PCA算法流程

4.3 案例演示

4.3.1 手工计算实例

4.3.2 使用 Python 实现

4.3.3 Open3D 计算点云法向量并显示

5.SVD分解法

6.知识点补充：已知平面方程，如何求法向量?

网站公告

今日签到

热门文章

最新发布