【计算方法与科学建模】矩阵特征值与特征向量的计算（二）：Jacobi 过关法及其Python实现（Jacobi 旋转法的改进）

  矩阵的特征值（eigenvalue）和特征向量（eigenvector）在很多应用中都具有重要的数学和物理意义。Jacobi 旋转法是一种用于计算对称矩阵特征值和特征向量的迭代方法，Jacobi 过关法是 Jacobi 旋转法的一种改进版本，其主要目的是减少计算工作和提高运行速度。

  本文将详细介绍Jacobi 过关法的基本原理和步骤，并给出其Python实现。

一、Jacobi 旋转法

  Jacobi 旋转法的每一次迭代中，需要选择一个非对角元素最大的位置，然后构造相应的旋转矩阵，进行相似变换，使得矩阵逐渐对角化。

• 对称矩阵是一个实数矩阵，其转置与自身相等。
• 对于一个方阵$A$，如果存在标量 $\lambda$ 和非零向量 $v$，使得 $Av=\lambda v$，那么 $\lambda$ 就是 $A$ 的特征值，$v$ 就是对应于 $\lambda$ 的特征向量。

1. 基本思想

  Jacobi 旋转法的基本思想是通过一系列的相似变换，逐步将对称矩阵对角化，使得非对角元素趋于零。这个过程中，特征值逐渐浮现在对角线上，而相应的特征向量也被逐步找到。下面是 Jacobi 旋转法的基本步骤：

1. 选择旋转角度： 选择一个旋转角度 θ，通常使得旋转矩阵中的非对角元素为零，从而实现对角化，通常选择非对角元素中绝对值最大的那个作为旋转的目标。

2. 构造旋转矩阵： 构造一个旋转矩阵 J，该矩阵为单位矩阵，只有对应于选择的非对角元素的位置上有两个非零元素，其余位置上为零。这两个非零元素的值由旋转角度 θ 决定，例如，对于 2x2 矩阵，旋转矩阵可以表示为：
   J = [ cos ⁡ ( θ ) − sin ⁡ ( θ ) sin ⁡ ( θ ) cos ⁡ ( θ ) ] J =

   [cos⁡(θ)−sin⁡(θ)sin⁡(θ)cos⁡(θ)]" role="presentation" style="position: relative;">[cos(θ)sin(θ)−sin(θ)cos(θ)]$$\left[\begin{array}{cc}\cos (\theta )& -\sin (\theta )\\ \sin (\theta )& \cos (\theta )\end{array}\right]$$
   J=[cos(θ)sin(θ)​−sin(θ)cos(θ)​]

3. 相似变换： 计算相似变换矩阵 $P$，即 $P^{T}AP$，其中 $A$ 是原始矩阵，$P$ 是旋转矩阵，计算过程如下：

P T A P = [ cos ⁡ ( θ ) sin ⁡ ( θ ) − sin ⁡ ( θ ) cos ⁡ ( θ ) ] T [ a 11 a 12 a 12 a 22 ] [ cos ⁡ ( θ ) − sin ⁡ ( θ ) sin ⁡ ( θ ) cos ⁡ ( θ ) ] P^TAP =

^T PTAP=[cos(θ)−sin(θ)​sin(θ)cos(θ)​]T[a11​a12​​a12​a22​​][cos(θ)sin(θ)​−sin(θ)cos(θ)​]

  通过矩阵相乘计算，我们可以得到 $P^{T}AP$ 中的非对角元素，假设这两个元素分别位于矩阵的 (1,2) 和 (2,1) 的位置。令 $a_{ij}$ 为这两个元素，即 $a_{ij}=a_{12}=a_{21}$。

  接下来，我们希望通过选择合适的 $\theta$使得 $a_{ij}$ 变为零，从而达到对角化的目的，即$a_{12}=a_{21}$，进一步可推导出

$$\theta =\frac{1}{2}\arctan \left(\frac{2\cdot a_{ij}}{a_{ii}-a_{jj}}\right)$$

• 若$a_{ii}=a_{jj}$，则使用$arccot$形式

4. 迭代： 重复步骤 1-3，直到矩阵 A 的非对角元素都趋于零或满足一定的精度要求。

5. 提取特征值和特征向量： 对角线上的元素即为矩阵 A 的特征值，而 P 中的列向量即为对应于这些特征值的特征向量。

2. 注意事项

  Jacobi 旋转法的优点是可以用于任意大小的对称矩阵，但其缺点是迭代次数较多，计算量较大。在实际应用中，通常会结合其他方法来提高计算效率。

二、Jacobi 过关法

  Jacobi 过关法（Jacobi’s threshold method）是 Jacobi 旋转法的一种改进版本，其主要目的是减少计算工作和提高运行速度。该方法通过动态调整阈值，并根据阈值对非对角元素进行选择性的旋转变换，以逐步对角化对称矩阵。

1. 基本思想

1. 计算非对角元素平方和： 对于对称矩阵 $A$，计算其非对角元素平方和，表示为 $u(A)$。然后取平方根，得到 $r(A)=\sqrt{u(A)}$。

2. 设定初始阈值 $\theta$： 预先设定一个初始阈值 ${\theta }_0$。

3. 扫描非对角元素： 对于 $a_{ij}$ 其中 $i\ne j$，扫描矩阵的上三角或下三角部分。

4. 进行选择性旋转变换： 对于绝对值大于当前阈值 $\theta $的非对角元素 $a_{ij}$，进行 Jacobi 旋转变换，将其旋转为零。旋转变换的具体步骤如下：

   • 选择旋转角度 $\varphi$，通常通过 $\tan (2\varphi )=\frac{2a_{ij}}{a_{ii}-a_{jj}}$计算。
   • 构造旋转矩阵 $J$，其中除了 $J_{ii}$ 和 $J_{jj}$ 外的元素都为零，而 $J_{ij}$ 和 $J_{ji}$ 元素由 $\cos (\varphi )$ 和 $-\sin (\varphi )$决定。
   • 执行相似变换 $A=J^{T}AJ$。

5. 调整阈值 $\theta$： 当所有非对角元素的绝对值都小于当前阈值 $\theta$ 时，缩小阈值，即 ${\theta }_{i+1}=\gamma \cdot {\theta }_i$，其中 $\gamma$ 是一个缩小因子。

6. 重复步骤 3-5： 重复上述步骤，直到满足某个收敛条件，例如 ${\theta }_k<ϵ$，其中 $ϵ$是一个很小的正数。

2. 注意事项

  通过不断调整阈值并选择性地进行旋转变换，Jacobi 过关法逐渐减小非对角元素的绝对值，以达到更好的数值稳定性。这种方法的优点在于，通过智能地选择非对角元素进行变换，可以有效减少迭代次数，提高计算效率。

三、Python实现

import numpy as np


def jacobi_threshold_method(A, epsilon=1e-10, gamma=0.9):
    n = A.shape[0]
    theta = np.sqrt(np.sum(np.abs(np.triu(A, k=1)) ** 2))
    eigenvectors = np.eye(n)

    while theta > epsilon:
        for i in range(n):
            for j in range(i + 1, n):
                if np.abs(A[i, j]) > theta:
                    # 计算旋转角度
                    phi = 0.5 * np.arctan2(2 * A[i, j], A[i, i] - A[j, j])

                    # 构造旋转矩阵
                    J = np.eye(n)
                    J[i, i] = J[j, j] = np.cos(phi)
                    J[i, j] = -np.sin(phi)
                    J[j, i] = np.sin(phi)

                    # 执行相似变换
                    A = np.dot(np.dot(J.T, A), J)

                    # 更新特征向量
                    eigenvectors = np.dot(eigenvectors, J)

        # 缩小阈值
        theta *= gamma

    # 提取特征值和特征向量
    eigenvalues = np.diag(A)

    return eigenvalues, eigenvectors


# 示例矩阵
A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])

# 执行 Jacobi 过关法
eigenvalues, eigenvectors = jacobi_threshold_method(A)

print("特征值:", eigenvalues)
print("特征向量:")
np.set_printoptions(precision=4, suppress=True)
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

迭代过程（调试）

• 第一次：

………

• final：