Jacobi迭代法求对称矩阵特征值及特征向量

我将为您讲解Jacobi迭代法求解对称矩阵特征值和特征向量的完整过程。这个经典算法通过一系列平面旋转逐步将对称矩阵对角化。

一、问题描述

给定一个 \(n \times n\) 的实对称矩阵 \(A\)，我们需要：

1. 计算所有特征值 \(\lambda_1, \lambda_2, \dots, \lambda_n\)
2. 计算对应的特征向量 \(v_1, v_2, \dots, v_n\)

关键特点：

• 对称矩阵特征值都是实数
• 特征向量可以构成正交基
• Jacobi方法通过迭代产生正交矩阵序列，逐步将A对角化

二、核心思想

基本思路

寻找一系列正交矩阵 \(J_k\)（平面旋转矩阵），使得：

\[A_{k+1} = J_k^T A_k J_k \]

其中 \(A_0 = A\)。当 \(k \to \infty\) 时，\(A_k\) 趋于对角矩阵 \(D\)：

\[D = \lim_{k \to \infty} A_k = \text{diag}(\lambda_1, \dots, \lambda_n) \]

累积的正交变换：

\[V = J_1 J_2 \cdots J_k \cdots \]

的列就是特征向量。

为什么选择平面旋转？

平面旋转矩阵 \(J(p, q, \theta)\) 在 \((p, q)\) 平面旋转角度 \(\theta\)，其他坐标不变：

\[J_{pp} = J_{qq} = \cos\theta,\quad J_{pq} = \sin\theta,\quad J_{qp} = -\sin\theta \]

其他对角元为1，非对角元为0。这种变换保持正交性且只改变特定行和列。

三、单次Jacobi旋转

1. 选择旋转元素

对于当前矩阵 \(A^{(k)}\)，选择非对角元中绝对值最大的：

\[|a_{pq}^{(k)}| = \max_{i

为什么选最大的？因为这样能使非对角元平方和减小最快。

2. 计算旋转角 θ

为了消去 \(a_{pq}\)，需要解：

\[\tan 2\theta = \frac{2a_{pq}}{a_{pp} - a_{qq}} \]

实际计算中为避免溢出，我们使用更稳定的公式：

设 \(\tau = \frac{a_{qq} - a_{pp}}{2a_{pq}}\)
则 \(t = \tan\theta\) 满足 \(t^2 + 2\tau t - 1 = 0\)
取绝对值较小的根：

\[t = \frac{\text{sign}(\tau)}{|\tau| + \sqrt{1+\tau^2}} \]

然后：

\[\cos\theta = \frac{1}{\sqrt{1+t^2}},\quad \sin\theta = t\cos\theta \]

3. 应用旋转

旋转矩阵 \(J = J(p, q, \theta)\) 作用于 \(A^{(k)}\)：

\[A^{(k+1)} = J^T A^{(k)} J \]

仅修改第 p, q 行和列（利用对称性）：

• \(a_{pp}^{(k+1)} = a_{pp}^{(k)} - t \cdot a_{pq}^{(k)}\)
• \(a_{qq}^{(k+1)} = a_{qq}^{(k)} + t \cdot a_{pq}^{(k)}\)
• \(a_{pq}^{(k+1)} = a_{qp}^{(k+1)} = 0\)

对于 \(i \neq p, q\)：

\[\begin{aligned} a_{ip}^{(k+1)} &= a_{pi}^{(k+1)} = a_{ip}^{(k)}\cos\theta - a_{iq}^{(k)}\sin\theta \\ a_{iq}^{(k+1)} &= a_{qi}^{(k+1)} = a_{ip}^{(k)}\sin\theta + a_{iq}^{(k)}\cos\theta \end{aligned} \]

4. 更新特征向量矩阵

累积变换：

\[V^{(k+1)} = V^{(k)} J \]

其中 \(V^{(0)} = I\)。每次只需更新第 p, q 列：

\[\begin{aligned} v_{ip}^{(k+1)} &= v_{ip}^{(k)}\cos\theta - v_{iq}^{(k)}\sin\theta \\ v_{iq}^{(k+1)} &= v_{ip}^{(k)}\sin\theta + v_{iq}^{(k)}\cos\theta \end{aligned} \]

四、完整算法流程

算法步骤

1. 初始化：

   • \(A^{(0)} = A\)
   • \(V^{(0)} = I_n\)
   • 设置容差 \(\epsilon\)（如 \(10^{-10}\)）

2. 计算非对角元范数：

\[ \text{off}(A) = \sqrt{\sum_{i

3. 迭代直到收敛（\(\text{off}(A) < \epsilon\)）：
   a) 找出最大非对角元 \(|a_{pq}|\)
   b) 计算旋转角 \(\theta\)
   c) 更新矩阵 \(A\)
   d) 更新特征向量矩阵 \(V\)
   e) 重新计算或更新 off-norm

4. 输出结果：

   • 特征值：\(\lambda_i = a_{ii}\)（最终的对角元）
   • 特征向量：\(V\) 的列

五、重要细节与优化

1. 循环Jacobi与经典Jacobi

• 经典Jacobi：每次选最大非对角元（需要 \(O(n^2)\) 搜索）
• 循环Jacobi：按固定顺序遍历所有 \((p, q)\) 对
  • 行循环：\((1,2), (1,3), \dots, (1,n), (2,3), \dots\)
  • 一次完整遍历称为一次"sweep"

2. 收敛性分析

设 \(S = \sum_{i\neq j} a_{ij}^2\) 为非对角元平方和。
一次旋转后：

\[S^{(k+1)} = S^{(k)} - 2a_{pq}^2 \]

因为每次消去最大的 \(a_{pq}\)，有 \(a_{pq}^2 \geq \frac{2S^{(k)}}{n(n-1)}\)，所以：

\[S^{(k+1)} \leq \left(1 - \frac{4}{n(n-1)}\right) S^{(k)} \]

线性收敛，但实际超线性收敛。

3. 数值稳定性

• 旋转公式设计避免小分母
• 使用双精度避免累积误差
• 特征向量通过累积旋转得到，而非最后求解

4. 阈值Jacobi

为减少小元素的无效旋转：
仅当 \(|a_{pq}| \geq \text{threshold}\) 时才旋转。
阈值可逐渐减小，如：

\[\text{threshold} = \frac{\text{off}(A)}{n\sqrt{2}} \]

六、算法伪代码

def jacobi_eigen(A, max_iter=1000, tol=1e-10):
    n = A.shape[0]
    V = np.eye(n)
    off_norm = np.sqrt(np.sum(np.triu(A, 1)**2) * 2)
    
    for iteration in range(max_iter):
        if off_norm < tol:
            break
            
        # 寻找最大非对角元 (p, q)
        p, q = np.unravel_index(np.argmax(np.abs(np.triu(A, 1))), A.shape)
        
        # 计算旋转角
        if A[p, p] == A[q, q]:
            theta = np.pi/4
        else:
            tau = (A[q, q] - A[p, p]) / (2 * A[p, q])
            t = np.sign(tau) / (abs(tau) + np.sqrt(1 + tau**2))
            c = 1 / np.sqrt(1 + t**2)
            s = t * c
        
        # 更新A
        app = A[p, p]
        aqq = A[q, q]
        apq = A[p, q]
        
        A[p, p] = app - t * apq
        A[q, q] = aqq + t * apq
        A[p, q] = A[q, p] = 0
        
        for i in range(n):
            if i != p and i != q:
                aip = A[i, p]
                aiq = A[i, q]
                A[i, p] = A[p, i] = c * aip - s * aiq
                A[i, q] = A[q, i] = s * aip + c * aiq
        
        # 更新特征向量
        for i in range(n):
            vip = V[i, p]
            viq = V[i, q]
            V[i, p] = c * vip - s * viq
            V[i, q] = s * vip + c * viq
        
        # 更新off-norm（或重新计算）
        off_norm = np.sqrt(off_norm**2 - 2*apq**2)
    
    eigenvalues = np.diag(A)
    eigenvectors = V
    return eigenvalues, eigenvectors

七、复杂度与特点

时间复杂度

• 每次旋转：\(O(n)\) 操作
• 一次sweep（所有对）：\(O(n^3)\)
• 总迭代：通常需要 \(O(\log(1/\epsilon))\) 次sweep
• 总体：\(O(n^3 \log(1/\epsilon))\)

空间复杂度

• \(O(n^2)\) 存储矩阵和特征向量

优点

1. 简单易懂，实现直接
2. 数值稳定（只有正交变换）
3. 能同时得到所有特征值和特征向量
4. 精度高（特征向量正交性好）

缺点

1. 对大型矩阵较慢（相比Lanczos/QR）
2. 不适合稀疏矩阵（会破坏稀疏性）
3. 收敛可能较慢，尤其特征值相近时

八、应用场景

1. 中小型稠密对称矩阵（n < 1000）
2. 需要高精度特征向量的场合
3. 并行计算（可同时处理多个旋转对）
4. 教学和理解特征值问题的经典案例

总结

Jacobi方法通过一系列精心选择的平面旋转，逐步将对称矩阵对角化。虽然对于大规模问题效率不如现代算法，但其优雅的几何直观性、数值稳定性和同时提供特征向量的能力，使其在中小规模问题中仍有重要价值。理解Jacobi方法是掌握矩阵特征值计算的基础。