当前位置：首页 > news >正文

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

news 来源：原创 2025/5/18 6:51:40

参考笔记：

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

1.前言

2.线性回归法

2.1 模型假设

2.2 定义误差函数

2.3 求偏导并解方程

2.4 解方程

2.5 案例演示

2.5.1 手工计算实例

2.5.2 使用Python实现

3.梯度下降法

4.PCA法

5.SVD分解法

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~1label=f'拟合平面')plt.show()

运行结果：

3.梯度下降法

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

4.PCA法拟合平面

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

另外，利用PCA法拟合平面还可以求出所拟合平面的法向量

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)$ ，将其拼成一个矩阵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}$

② 中心化数据

求出x维的均值 $\color{red}\bar{x}$ ，y维的均值 $\color{red}\bar{y}$ ，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.2 案例演示

4.2.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.2.2 使用 Python 实现

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

import numpy as np
import matplotlib.pyplot as pltplt.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~1label=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()

运行结果：