唐承乾

Jacobi迭代方法概述

先回顾下什么是特征值与特征向量：

对于一个$n \times n$ 的方阵 $A$，如果存在一个非零向量 $\mathbf{x}$ 和一个标量 $\lambda$，使得

则称 $\lambda$为矩阵 $A$ 的特征值，$\mathbf{x}$ 称为对应的特征向量。 如果 $A$ 是对称矩阵（即 $A = A^T$），那么它具有以下重要性质：

• 所有特征值都是实数。

• 不同特征值对应的特征向量彼此正交。

• 存在一个正交矩阵 $P$ 将 $A$对角化，即

  $$P^T A P = D$$

  其中 $D$ 是对角矩阵，其对角元素即为 $A$ 的特征值。

Jacobi迭代方法是一种用于对称矩阵特征值分解的数值方法。其核心思想是通过一系列旋转矩阵，逐步将原矩阵对角化。每次迭代选取矩阵中最大的非对角元素，构造一个旋转矩阵，使得该元素在旋转后变为零。重复这一过程，直到矩阵近似为对角矩阵，此时对角线上的元素即为特征值，旋转矩阵的乘积给出特征向量。领悟不到没关系，细看下算法实现步骤和两个例子就可以明白了。

在实数域上，对于对称矩阵 A，Jacobi方法通常能够稳定收敛。当矩阵的最大非对角元素小到一定阈值时，便可认为矩阵已基本对角化，停止迭代。

算法实现

1. 初始化：设定对称矩阵 $A$ 和单位矩阵 $V$ 作为特征向量的初始矩阵。
2. 查找最大非对角元素：在矩阵 $A$ 中找到绝对值最大的非对角元素 $A_{pq}$。
3. 计算旋转角度： $$\theta = \frac{1}{2} \arctan\left(\frac{2A_{pq}}{A_{qq} - A_{pp}}\right)$$ 计算旋转矩阵的正弦和余弦值： $$c = \cos(\theta), \quad s = \sin(\theta)$$
4. 构造旋转矩阵：生成一个单位旋转矩阵 $J$，其中仅在 $(p,p)$，$(p,q)$，$(q,p)$，$(q,q)$ 位置有非单位元素： $$J = \begin{pmatrix} \ddots & & & & \\ & c & \cdots & s & \\ & \vdots & 1 & \vdots & \\ & -s & \cdots & c & \\ & & & & \ddots \end{pmatrix}$$
5. 更新矩阵和特征向量： $$A' = J^T A J, \quad V' = V J$$
6. 检查停止条件：如果所有非对角元素的绝对值均小于预设阈值 $\epsilon$，则停止迭代；否则，返回步骤 2。

2x2矩阵的例子

考虑矩阵 $A = \begin{pmatrix} 4 & 1 \\ 1 & 3 \end{pmatrix}$。

1. 初始化： $$V = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$
2. 查找最大非对角元素：$A_{12} = 1$，唯一的非对角元素。
3. 计算旋转角度： $$\theta = \frac{1}{2} \arctan\left(\frac{2 \times 1}{3 - 4}\right) = \frac{1}{2} \arctan(-2) \approx -0.4636 \text{ 弧度}$$ $$c = \cos(-0.4636) \approx 0.8944, \quad s = \sin(-0.4636) \approx -0.4472$$
4. 构造旋转矩阵： $$J = \begin{pmatrix} 0.8944 & -0.4472 \\ 0.4472 & 0.8944 \end{pmatrix}$$
5. 更新矩阵： $$A' = J^T A J = \begin{pmatrix} 4.2361 & 0 \\ 0 & 2.7639 \end{pmatrix}$$ 更新特征向量： $$V' = V J = \begin{pmatrix} 0.8944 & -0.4472 \\ 0.4472 & 0.8944 \end{pmatrix}$$
6. 检查停止条件：所有非对角元素为 0，满足条件，迭代结束。

特征值为 4.2361 和 2.7639，特征向量分别为 $\begin{pmatrix} 0.8944 \\ 0.4472 \end{pmatrix}$ 和 $\begin{pmatrix} -0.4472 \\ 0.8944 \end{pmatrix}$。

3x3矩阵的例子

考虑矩阵 $A = \begin{pmatrix} 4 & 1 & 2 \\ 1 & 3 & 1 \\ 2 & 1 & 3 \end{pmatrix}$。

1. 初始化： $$V = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}$$
2. 查找最大非对角元素：$A_{13} = 2$ 是最大的非对角元素。
3. 计算旋转角度： $$\theta = \frac{1}{2} \arctan\left(\frac{2 \times 2}{3 - 4}\right) = \frac{1}{2} \arctan(-4) \approx -1.1903 \text{ 弧度}$$ $$c = \cos(-1.1903) \approx 0.3624, \quad s = \sin(-1.1903) \approx -0.9320$$
4. 构造旋转矩阵（对第1和第3列进行旋转）： $$J = \begin{pmatrix} 0.3624 & 0 & -0.9320 \\ 0 & 1 & 0 \\ 0.9320 & 0 & 0.3624 \end{pmatrix}$$
5. 更新矩阵： $$A' = J^T A J \approx \begin{pmatrix} 5.0 & 1.7321 & 0 \\ 1.7321 & 3 & 1 \\ 0 & 1 & 1 \end{pmatrix}$$ 更新特征向量： $$V' = V J = \begin{pmatrix} 0.3624 & 0 & -0.9320 \\ 0 & 1 & 0 \\ 0.9320 & 0 & 0.3624 \end{pmatrix}$$
6. 下一步迭代：重复上述步骤，继续查找更新后的最大非对角元素，直至所有非对角元素小于阈值。

Jacobi迭代方法在对称矩阵下保证收敛，即逐步将矩阵对角化。停止条件通常设定为所有非对角元素的绝对值均小于预设的阈值 $\epsilon$，如 $10^{-6}$。收敛速度依赖于矩阵的初始形式和选择旋转顺序，但在实践中通常能较快达到所需精度。

MATLAB代码示例

以下是使用MATLAB实现Jacobi迭代方法的简洁代码。

1
function [D, V] = jacobi_eigen(A, tol, max_iter)
2
    % 检查是否为方阵和对称矩阵
3
    [n, m] = size(A);
4
    V = eye(n);
5
    iter = 0;
6
    while true
7
        iter = iter + 1;
8
        % 查找最大的非对角元素
9
        off = A - diag(diag(A));
10
        [max_val, ind] = max(abs(off(:)));
11
        [p, q] = ind2sub(size(A), ind);
12
        if max_val < tol || iter > max_iter
13
            break;
14
        end
15
        if A(p,p) == A(q,q)
16
            theta = pi/4;
17
        else
18
            theta = 0.5 * atan(2*A(p,q)/(A(q,q) - A(p,p)));
19
        end
20
        c = cos(theta);
21
        s = sin(theta);
22
        % 构造旋转矩阵
23
        J = eye(n);
24
        J(p,p) = c;
25
        J(q,q) = c;
26
        J(p,q) = s;
27
        J(q,p) = -s;
28
        % 更新A和V
29
        A = J' * A * J;
30
        V = V * J;
31
    end
32
    D = A;
33
end
34

35
% 2x2矩阵示例
36
A2 = [4, 1; 1, 3];
37
tol = 1e-6;
38
max_iter = 100;
39
[D2, V2] = jacobi_eigen(A2, tol, max_iter);
40
disp('2x2矩阵特征值:');
41
disp(diag(D2));
42
disp('2x2矩阵特征向量:');
43
disp(V2);
44

45
% 3x3矩阵示例
46
A3 = [4, 1, 2; 1, 3, 1; 2, 1, 3];
47
[D3, V3] = jacobi_eigen(A3, tol, max_iter);
48
disp('3x3矩阵特征值:');
49
disp(diag(D3));
50
disp('3x3矩阵特征向量:');
51
disp(V3);

matlab计算结果如下