拟合平面再思考

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 - centroid

    xx = 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 * yz
    det_y = xx * zz - xz * xz
    det_z = xx * yy - xy * xy

    if 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

3。特征值求解原理

3。1 jacobi

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

1. 选择旋转角度： 在每一次迭代中，Jacobi 方法会选择一个最大的非对角线元素，然后计算旋转角度。旋转角度的选择有多种方法，一种常见的是选取使得旋转后的非对角线元素为零的角度。

2. 构造旋转矩阵： 利用选择的旋转角度，构造一个旋转矩阵。这个旋转矩阵是一个正交矩阵，它的转置等于它的逆。

3. 相似变换： 将原矩阵通过相似变换，即左乘和右乘旋转矩阵，得到一个新的矩阵。这个新矩阵在对角线上的元素更接近特征值。

4. 迭代： 重复以上步骤，直到矩阵的非对角线元素足够小，或者达到预定的迭代次数。

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

3。2对角化

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

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

[ P^{-1}AP = D ]

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

对角化的步骤如下：

1. 计算特征值： 找到矩阵 (A) 的特征值 (\lambda).

2. 计算特征向量： 对于每个特征值，找到对应的特征向量 (\mathbf{v}).

3. 构造矩阵 (P)： 将特征向量按列排列成矩阵 (P).

4. 对角化： 计算 (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}) 是对角线上的元素。

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

1. 幂次运算简化： 对角矩阵的幂次运算非常简单。如果 (D) 是对角矩阵，那么 (D^k) 的每个对角元素就是 (d_{i}^k)。

2. 特征值和特征向量： 对角矩阵的特征值就是它的对角元素，而每个对角元素本身就是对应的特征向量。

3. 可逆性： 对角矩阵是可逆的当且仅当每个对角元素都不为零。

4. 乘法和加法： 对角矩阵的乘法和加法操作都很简单，因为大部分元素都是零。

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

3。4 python实现

import numpy as np

def 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 = 0
    while 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 / 4
        else:
            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_matrix
        eigenvectors = eigenvectors @ rotation_matrix

        iter_count += 1

    return 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)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64