Jacobi迭代方法矩阵特征值分解

Jacobi迭代方法矩阵特征值分解

Jacobi迭代方法概述

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

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

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

• 所有特征值都是实数。

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

• 存在一个正交矩阵 $P$ 将 $A$对角化，即
  $$P^{T}AP=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 ),s=\sin (\theta )$$

4. 构造旋转矩阵：生成一个单位旋转矩阵 $J$，其中仅在 $(p,p)$，$(p,q)$，$(q,p)$，$(q,q)$ 位置有非单位元素：
   $$J=\left(\begin{array}{ccccc}\ddots & & & & \\ & c& \cdots & s& \\ & ⋮& 1& ⋮& \\ & -s& \cdots & c& \\ & & & & \ddots \end{array}\right)$$

5. 更新矩阵和特征向量：
   $$A^{\prime }=J^{T}AJ,V^{\prime }=VJ$$

6. 检查停止条件：如果所有非对角元素的绝对值均小于预设阈值 $ϵ$，则停止迭代；否则，返回步骤 2。

2x2矩阵的例子

考虑矩阵 $A=\left(\begin{array}{cc}4& 1\\ 1& 3\end{array}\right)$。

1. 初始化：
   $$V=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$

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,s=\sin (-0.4636)\approx -0.4472$$

4. 构造旋转矩阵：
   $$J=\left(\begin{array}{cc}0.8944& -0.4472\\ 0.4472& 0.8944\end{array}\right)$$

5. 更新矩阵：
   $$A^{\prime }=J^{T}AJ=\left(\begin{array}{cc}4.2361& 0\\ 0& 2.7639\end{array}\right)$$
   更新特征向量：
   $$V^{\prime }=VJ=\left(\begin{array}{cc}0.8944& -0.4472\\ 0.4472& 0.8944\end{array}\right)$$

6. 检查停止条件：所有非对角元素为 0，满足条件，迭代结束。

特征值为 4.2361 和 2.7639，特征向量分别为
$$\left(\begin{array}{c}0.8944\\ 0.4472\end{array}\right)$$
和
$$\left(\begin{array}{c}-0.4472\\ 0.8944\end{array}\right)$$。

3x3矩阵的例子

考虑矩阵 $A=\left(\begin{array}{ccc}4& 1& 2\\ 1& 3& 1\\ 2& 1& 3\end{array}\right)$。

1. 初始化：
   $$V=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)$$

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,s=\sin (-1.1903)\approx -0.9320$$

4. 构造旋转矩阵（对第1和第3列进行旋转）：
   $$J=\left(\begin{array}{ccc}0.3624& 0& -0.9320\\ 0& 1& 0\\ 0.9320& 0& 0.3624\end{array}\right)$$

5. 更新矩阵：
   $$A^{\prime }=J^{T}AJ\approx \left(\begin{array}{ccc}5.0& 1.7321& 0\\ 1.7321& 3& 1\\ 0& 1& 1\end{array}\right)$$
   更新特征向量：
   $$V^{\prime }=VJ=\left(\begin{array}{ccc}0.3624& 0& -0.9320\\ 0& 1& 0\\ 0.9320& 0& 0.3624\end{array}\right)$$

6. 下一步迭代：重复上述步骤，继续查找更新后的最大非对角元素，直至所有非对角元素小于阈值。

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

MATLAB代码示例

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

function [D, V] = jacobi_eigen(A, tol, max_iter)
    % 检查是否为方阵和对称矩阵
    [n, m] = size(A);
    V = eye(n);
    iter = 0;
    while true
        iter = iter + 1;
        % 查找最大的非对角元素
        off = A - diag(diag(A));
        [max_val, ind] = max(abs(off(:)));
        [p, q] = ind2sub(size(A), ind);
        if max_val < tol || iter > max_iter
            break;
        end
        if A(p,p) == A(q,q)
            theta = pi/4;
        else
            theta = 0.5 * atan(2*A(p,q)/(A(q,q) - A(p,p)));
        end
        c = cos(theta);
        s = sin(theta);
        % 构造旋转矩阵
        J = eye(n);
        J(p,p) = c;
        J(q,q) = c;
        J(p,q) = s;
        J(q,p) = -s;
        % 更新A和V
        A = J' * A * J;
        V = V * J;
    end
    D = A;
end

% 2x2矩阵示例
A2 = [4, 1; 1, 3];
tol = 1e-6;
max_iter = 100;
[D2, V2] = jacobi_eigen(A2, tol, max_iter);
disp('2x2矩阵特征值:');
disp(diag(D2));
disp('2x2矩阵特征向量:');
disp(V2);

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

matlab计算结果如下