一、幂法和反幂法

1、幂法

由于求解特征方程的计算量巨大，所以在求解特征值时，要从特征值的定义

入手。设 阶实矩阵 的特征值 满足：

且与 相应的特征向量 线性无关。给定初始向量 ，由迭代公式：

产生向量序列 ，当 充分大时，有 ，相应的特征向量为 。

证明当充分大时，上述结论成立：

设 ，且 线性无关，故必存在 个不全为零的数 ，使得：

得：

即：

设 , 由 得：

于是：

故只要 充分大, 就有：

因此, 可把 作为与 相应的特征向量的近似。由公式：

即可得到：

利用上述结论，可以不断迭代，得到：

注意到其中在较大的情况下会出现两种状况：

1. 如果的范数大于，则会导致计算过程中溢出。
2. 如果的范数小于，则会导致在计算过程中变成。

根据上面的推导过程，以及矩阵2-范数的定义，为了避免这种情形，实际计算时每次迭代所求得的向量都要归一化，因此实际计算使用公式为：

当我们计算的结果小于判定值时，结束循环，判定为其已经收敛，得到特征向量，在此之后，一般不做除法来获得特征值，而是使用瑞利商：

这样所得的结果更加稳定。

2、反幂法

在实际应用是，我们需要寻找一个最接近于数的特征值，可以通过求解的最小特征值（绝对值最小），从而将所得结果加后获得答案。

设 为 阶非奇异矩阵， 为 的特征值与相应的特征向量，则 的特征值是 的特征值的倒数，而相应的特征向量不变，即

因此，若对矩阵 用幂法，即可计算出 的按模最大的特征值，其倒数恰为 的按模最小的特征值。这就是反幂法的基本思想。

由于的计算比较麻烦，而且往往不能保持矩阵的一些好的性质（如稀疏性），因此反幂法在实际运算时以求解方程组：

代替，同时在每次迭代的过程中做归一化处理。

3、代码

下面是使用反幂法求解特定特征值（距离给出 lambda0 最近的特征值）的MATLAB代码：

1
function [lambda, v] = solveEigenpair(A, lambda0)
2
    n = size(A, 1);
3
    B = A - lambda0 * eye(n);
4
    [L, U, P] = lu(B);
5
    v = rand(n, 1);
6
    v = v / norm(v);
7
    tol = 1e-12;
8
    while true
9
        v_new = U \ (L \ (P * v));
10
        v_new = v_new / norm(v_new);
11
        lambda = v_new' * A * v_new;
12
        if norm(A * v_new - lambda * v_new) < tol
13
            break;
14
        end
15
        v = v_new;
16
    end
17
end

二、雅可比（Jacobi）方法

Jacobi 方法是求实对称矩阵的全部特征值及相应特征向量的一种方法 。它基于以下两个数学结论：

1. 正交相似变换：任意实对称矩阵 均存在正交矩阵 ，使得 ，其中 为特征值， 的列向量为特征向量。
2. Frobenius 范数不变性：在正交变换下，矩阵元素的平方和保持不变 。设 ，则 。

该方法的核心思想是通过一系列正交变换（旋转变换），不断减小非对角元素的平方和，使矩阵趋向于对角阵。下面介绍旋转变换矩阵 ，旋转变换矩阵是将单位阵的第 行/列交点处的元素替换为三角函数值得到：

选取旋转角 使得变换后 。若 ，则取：

雅各比方法运算量大，且不能保持矩阵的特殊形状，所以该方法只适合求解中小型稠密矩阵。

三、QR 方法

1、基本 QR 算法

QR 算法是计算中小型矩阵全部特征值最有效的方法之一 。其基本思想是利用矩阵的 QR 分解进行相似变换迭代，可以证明，任意非奇异实矩阵都可以分解成一个正交矩阵和一个上三角矩阵的乘积，并且当的对角元符号确定时，分解时唯一的，其分解格式为：

其中 为正交阵， 为上三角阵 。迭代产生的序列 相似于原矩阵，并最终收敛于上三角形（或分块上三角形）矩阵。因为上三角阵的主对角元是该矩阵的特征值，故充分大时，的主对角元可以看作是原矩阵的特征值近似。

使用施密特（Schmidt）正交化方法对 阶非奇异实矩阵 进行 QR 分解的通式推导过程如下：

设 为 阶非奇异实矩阵，记为列向量形式：

其中 。利用施密特正交化构造正交向量组 ，对于第一列向量 ，直接进行单位化：

对于后续列向量 ，先求正交向量 ，从 中减去其在已求出的正交向量 上的投影：

得到长度为 1 的正交向量 ：

将上述迭代公式改写为 的线性组合形式：

根据该线性组合关系，矩阵 可以写成正交矩阵 与上三角矩阵 的乘积：

2、Householder 变换

为了减少运算量，通常先通过 Householder 变换将一般矩阵化为拟上三角形（Hessenberg 矩阵）。设 ，则称 为 Householder 矩阵。 是实对称的正交矩阵（），且 。它可以将任一非零向量 变换为单位向量的倍数（即实现坐标轴旋转/反射），用于矩阵消零。

该矩阵存在一个重要性质：设 为 中任意两个非零向量，且 ，则存在 Householder 矩阵 ，使得：

该定理表明，对于任一非零向量 ，都可以构造一个 Householder 变换，将其变成事先给定的单位向量的倍数。在数值计算中，通常取 （单位坐标向量），使 变换后只有一个分量不为零。

为了避免在 与 方向接近时产生较大的舍入误差，实际计算中总取以下通式构造 ：

3、化一般矩阵为拟上三角矩阵

拟上三角形矩阵（上海森堡阵）指主对角线下方的元素从第二行开始全为零的矩阵 。其形式如下：

首先构造第一个 Householder 矩阵 ，选取 矩阵第一列对角线以下的元素组成向量 。根据定理，构造一个 阶的 Householder 矩阵 ，使得变换后的向量 只有第一个分量不为零，即 。将 嵌入到 阶单位阵中：

计算 ，此时所得矩阵的第一列除前两个元素外，其余元素均为零。对于 ，依次构造 。 的形式为：

其中 是 阶单位阵， 是对当前矩阵第 列对角线以下元素进行处理的 Householder 矩阵。经过 次正交相似变换后，得到拟上三角形矩阵 ：

4、拟上三角矩阵的QR分解

由于拟上三角形矩阵 的特殊形状（仅在次对角线上有非零元），通常使用 个旋转变换（即 Givens 变换）将其化为上三角形矩阵，从而实现 QR 分解 。通过一系列旋转矩阵 依次消去拟上三角矩阵 的次对角线元素 。

取旋转矩阵 消去 ：

计算参数：，，。变换，此时 。

循环操作，假设前 步已完成，取 消去 的次对角线元 。计算参数，， 。变换 。经过最多 次旋转变换后，矩阵变为上三角矩阵 ：

由于旋转矩阵均为正交矩阵，根据正交阵的性质：

其中， 仍为正交矩阵。完成这一过程的运算量约为，比一般矩阵的QR分交少了一个数量级。

5、代码

下面为使用 Householder 矩阵优化后的QR分解的MATLAB代码：

1
function [Q, H] = hessenberg(A)
2
    n = size(A, 1);
3
    Q = eye(n);
4
    H = A;
5
    for i = 1:n-2
6
        x = H(i+1:n, i);
7
        e = zeros(n-i, 1);
8
        e(1) = 1;
9
        s = sign(x(1));
10
        if s == 0
11
            s = 1;
12
        end
13
        w = (x + norm(x, 2) * e * s) / norm(x + norm(x, 2) * e * s, 2);
14
        H0 = eye(n-i) - 2 * (w * w');
15
        H1 = eye(n);
16
        H1(i+1:n, i+1:n) = H0;
17
        Q = Q * H1;
18
        H = H1 * H * H1;
19
    end
20
end
21

22
function [Q, R] = getQR(H)
23
    n = size(H, 1);
24
    Q = eye(n);
25
    R = H;
26
    for i = 1:n-1
27
        a = R(i, i);
28
        b = R(i+1, i);
29
        if b == 0
30
            c = 1;
31
            s = 0;
32
        else
33
            r = hypot(a, b);
34
            c = a / r;
35
            s = b / r;
36
        end
37
        ri = R(i, i:n);
38
        rj = R(i+1, i:n);
39
        R(i, i:n)   = c * ri + s * rj;
40
        R(i+1, i:n) = -s * ri + c * rj;
41
        qi = Q(:, i);
42
        qj = Q(:, i+1);
43
        Q(:, i)   = c * qi + s * qj;
44
        Q(:, i+1) = -s * qi + c * qj;
45
    end
46
end
47

48
function values = eigen(A)
49
    n = size(A, 1);
50
    [~, H] = hessenberg(A);
51
    tol = 1e-5;
52
    while norm(H - diag(diag(H)), 'fro') >= tol
53
        mu = H(n, n);
54
        B = H(1:n, 1:n) - mu * eye(n);
55
        [Q, R] = getQR(B);
56
        H = R * Q + mu * eye(n);
57
    end
58
    values = real(diag(H));
59
end