拟合平面再思考

0。总结

之前思考了2种拟合平面的方式
1。ransan 随机取样拟合平面
2。特征值分解特征向量，最小特征值就是法向量

其中方法一误差较小，毕竟随机采样的方式可以忽略误差点。而特征值分解的方式如果误差较大，那么得到的结果非常不好

1。扩展

https://www.ilikebigbits.com/2015_03_04_plane_from_points.html
这篇文章说，通过求协方差矩阵就可以更好的拟合平面
这个文章的重点是他假设c=1 也就是c不能等于0。也就是平面的法向量的z轴不能为0
然后求得 a b d 这样就求得了平面方程
ax+by+cz+d=0

但是看代码发现他还要求的最大值？
这是为了避免假设不成立。所以他计算了3次。
假设a=1 ，假设b=1 ，假设c=1。通过这种方式求了3次。获取最大值保证法向量的xyz轴分量假设不为0的情况

然后求得法向量

但是这个法向量和通过协方差矩阵求的还是有差异。具体的取舍就看大家了

2。代码

def getplanefrompoints(points):centroid = np.mean(points, axis=0)r = points - centroidxx = np.sum(r[:, 0] * r[:, 0])xy = np.sum(r[:, 0] * r[:, 1])xz = np.sum(r[:, 0] * r[:, 2])yy = np.sum(r[:, 1] * r[:, 1])yz = np.sum(r[:, 1] * r[:, 2])zz = np.sum(r[:, 2] * r[:, 2])det_x = yy * zz - yz * yzdet_y = xx * zz - xz * xzdet_z = xx * yy - xy * xyif det_x > det_y and det_x > det_z:vecC = (det_x, xz * yz - xy * zz, xy * yz - xz * yy)elif det_y > det_z:vecC = (xz * yz - xy * zz, det_y, xy * xz - yz * xx)else:vecC = (xy * yz - xz * yy, xy * xz - yz * xx, det_z)vecC = vecC / np.linalg.norm(vecC)d = -vecC.dot(centroid)return vecC

3。特征值求解原理

3。1 jacobi

Jacobi 方法是一种迭代法，用于计算对称矩阵的特征值和特征向量。它的基本思想是通过一系列相似变换，逐步将对称矩阵对角化，使得非对角线上的元素逐渐趋于零。以下是 Jacobi 方法的基本原理：

选择旋转角度： 在每一次迭代中，Jacobi 方法会选择一个最大的非对角线元素，然后计算旋转角度。旋转角度的选择有多种方法，一种常见的是选取使得旋转后的非对角线元素为零的角度。
构造旋转矩阵： 利用选择的旋转角度，构造一个旋转矩阵。这个旋转矩阵是一个正交矩阵，它的转置等于它的逆。
相似变换： 将原矩阵通过相似变换，即左乘和右乘旋转矩阵，得到一个新的矩阵。这个新矩阵在对角线上的元素更接近特征值。
迭代： 重复以上步骤，直到矩阵的非对角线元素足够小，或者达到预定的迭代次数。

Jacobi 方法的优点是简单易懂，容易实现，并且适用于小型矩阵。然而，它的缺点是收敛速度相对较慢，尤其是对于大型矩阵。在实际应用中，更高效的算法（如QR分解、Lanczos 方法等）通常会被使用，以便更快地计算特征值和特征向量。

3。2对角化

对角化是线性代数中的一个重要概念，它涉及将一个矩阵转换为对角矩阵的过程。对角化使得矩阵的很多运算变得更加简单，特别是计算矩阵的幂次、矩阵指数和矩阵函数等操作。

对于一个 (n \times n) 的方阵 (A)，如果存在一个可逆矩阵 (P) 使得 (P^{-1}AP) 是一个对角矩阵，即：

[ P^{-1}AP = D ]

其中，(D) 是一个对角矩阵，其对角线上的元素是 (A) 的特征值。这样的矩阵 (P) 被称为 (A) 的特征向量矩阵，而 (D) 则是 (A) 的特征值矩阵。

对角化的步骤如下：

计算特征值： 找到矩阵 (A) 的特征值 (\lambda).
计算特征向量： 对于每个特征值，找到对应的特征向量 (\mathbf{v}).
构造矩阵 (P)： 将特征向量按列排列成矩阵 (P).
对角化： 计算 (P^{-1}AP).

对角化的好处在于，对角矩阵的幂次运算非常简单，因为对角矩阵的幂次就是每个对角元素的幂次。这对于某些数值计算和线性代数的问题是非常有用的。

在实际应用中，对角化不一定对所有矩阵都可行。可对角化的矩阵被称为可对角化矩阵，而对于某些矩阵，它们可能不可对角化。

3。3对角矩阵

对角矩阵是一种特殊形式的矩阵，其中除了主对角线上的元素之外，所有其他元素都为零。一个 (n \times n) 的对角矩阵可以表示为：

[
D = \begin{bmatrix}
d_{1} & 0 & \cdots & 0 \
0 & d_{2} & \cdots & 0 \
\vdots & \vdots & \ddots & \vdots \
0 & 0 & \cdots & d_{n}
\end{bmatrix}
]

其中 (d_{1}, d_{2}, \ldots, d_{n}) 是对角线上的元素。

对角矩阵具有一些有趣的性质，使得它们在线性代数中非常有用。以下是一些关于对角矩阵的重要性质：

幂次运算简化： 对角矩阵的幂次运算非常简单。如果 (D) 是对角矩阵，那么 (D^k) 的每个对角元素就是 (d_{i}^k)。
特征值和特征向量： 对角矩阵的特征值就是它的对角元素，而每个对角元素本身就是对应的特征向量。
可逆性： 对角矩阵是可逆的当且仅当每个对角元素都不为零。
乘法和加法： 对角矩阵的乘法和加法操作都很简单，因为大部分元素都是零。

对角矩阵在许多数学和工程应用中经常出现，特别是在对角化矩阵的概念中。对角化是将一个矩阵转换为对角矩阵的过程，便于处理矩阵的一些运算。

3。4 python实现

import numpy as npdef jacobi_eigenvalue(A, tol=1e-10, max_iter=1000):"""Jacobi 方法求解对称矩阵的特征值和特征向量Parameters:- A: 对称矩阵- tol: 迭代收敛的容忍度- max_iter: 最大迭代次数Returns:- eigenvalues: 特征值数组- eigenvectors: 特征向量矩阵，每一列是对应的特征向量"""n = A.shape[0]eigenvalues = np.diag(A).copy()  # 初始时取矩阵的对角线元素作为特征值的初始估计eigenvectors = np.eye(n)  # 初始化特征向量矩阵为单位矩阵iter_count = 0while iter_count < max_iter:# 找到最大的非对角元素max_off_diag = np.max(np.abs(np.triu(A, k=1)))if max_off_diag < tol:break  # 收敛条件# 找到最大非对角元素的位置indices = np.where(np.abs(A) == max_off_diag)i, j = indices[0][0], indices[1][0]# 计算旋转角度if A[i, i] == A[j, j]:theta = np.pi / 4else:theta = 0.5 * np.arctan(2 * max_off_diag / (A[i, i] - A[j, j]))# 构造旋转矩阵rotation_matrix = np.eye(n)rotation_matrix[i, i] = np.cos(theta)rotation_matrix[j, j] = np.cos(theta)rotation_matrix[i, j] = -np.sin(theta)rotation_matrix[j, i] = np.sin(theta)# 进行相似变换A = rotation_matrix.T @ A @ rotation_matrixeigenvectors = eigenvectors @ rotation_matrixiter_count += 1return eigenvalues, eigenvectors# 生成一个对称矩阵作为例子
A = np.array([[4.0, -2.0, 2.0],[-2.0, 2.0, -4.0],[2.0, -4.0, 11.0]])eigenvalues, eigenvectors = jacobi_eigenvalue(A)
print("Eigenvalues:")
print(eigenvalues)print("\nEigenvectors:")
print(eigenvectors)