1. 幂法

算法原理：从非零初始向量​出发，由构造向量序列，k充分大时，​近似为主特征值​的特征向量，相邻迭代向量分量比收敛到​。

算法步骤：

a. 给定初始向量、误差限ε、最大迭代次数it_max，令k=1；

b. 计算；

c. 若，停止；

d. 未达最大迭代次数则k=k+1，转 b。最终​为​，为对应特征向量。

MATLAB 程序：

function [m, u, k] = pow(A, ep, it_max)
if nargin < 3    it_max = 100;    end
if nargin < 2    ep = 1e-5;    end
n = length(A);    u = ones(n, 1);    k = 0;    m1 = 0;
while k <= it_max
    v = A * u;    [vmax, i] = max(abs(v));    m = v(i);    u = v / m;
    if abs(m - m1) < ep    break;    end
    m1 = m;    k = k + 1;
end

数值实验

例 1：给定矩阵

用幂法求A的主特征值及特征向量：

• 先求全部特征值：

  >> A = [5 4 1 1; 4 5 1 1; 1 1 4 2; 1 1 2 4];
  >> eig(A)
  

• 用幂法求主特征值：

  >> ep = 1e-5;    it_max = 100;    [m, u, k] = pow(A, ep, it_max)
  

• 原点平移法加速（取p=3）：

  >> B = A - 3 * eye(4);    [m, u, k] = pow(B, ep, it_max)
  

  

• 迭代次数减少，加速显著。

2. 反幂法

算法原理：利用A按模最小特征值是按模最大特征值，对用幂法，得主特征值​，从而得A按模最小特征值​；也可用原点平移法加速。

算法步骤：

a. 给定初始向量、误差限ε、最大迭代次数it_max，令k=1；

b. 对A做 LU 分解A=LU；

c. 解得；

d. 计算​；

e. 若，停止；

f. 未达最大迭代次数则k=k+1，转 b。最终​为​，为对应特征向量。

MATLAB 程序

function [m, u, k] = pow_inv(A, ep, it_max)
if nargin < 3    it_max = 100;    end
if nargin < 2    ep = 1e-10;    end
n = length(A);    u = ones(n, 1);    k = 1;    m1 = 0;
[L, U, P] = lu(A);    L_1 = inv(L);    U_1 = inv(U);
v = U_1 * u;    [vmax, i] = max(abs(v));    m = v(i);    u = v / m;
while k <= it_max
    y = L_1 * P * u;    v = U_1 * y;
    [vmax, i] = max(abs(v));    m = v(i);    u = v / m;
    if abs(m - m1) < ep    break;    end
    m1 = m;    k = k + 1;
end
m = 1 / m;

数值实验（例 2，对例 1 矩阵）：

• 直接求按模最小特征值：

  >> A = [5 4 1 1; 4 5 1 1; 1 1 4 2; 1 1 2 4];    ep = 1e-10;    it_max = 100;
  >> [m, u, k] = pow_inv(A, ep, it_max)
  

                                                                                                             特征值 1，特征向量，迭代 15 次。

• 求接近 2 的特征值（取p=1.8）：

  >> p = 1.8;    B = A - p * eye(4);    [m, u, k] = pow_inv(B, ep, it_max)
  

                                                                                                  特征值0.2+1.8=2，特征向量，迭代 11 次。

• 求接近 5 的特征值（取p=4.5）：

  >> p = 4.5;    B = A - p * eye(4);    [m, u, k] = pow_inv(B, ep, it_max)
  

                                                                                                特征值0.5+4.5=5，特征向量，迭代 16 次。

矩阵A特征值 1,2,5,10，原点平移法中特征值分离越好，迭代越快。

3. 用正交相似变换约化一般矩阵为上黑森伯格矩阵

算法原理：对，存在初等反射矩阵​，使（上黑森伯格矩阵）。

算法步骤：

a. 输入方阵A；

b. 对k=1,2,⋯,n−2，依次做豪斯霍尔德反射变换（相似变换），将A约化为上黑森伯格矩阵。

MATLAB 程序：

function [Q, R] = hessenberg(A)
% 用householder正交相似变换将一般矩阵约化为上黑森伯格矩阵
[m, n] = size(A);
if m ~= n    error('A非方阵');    end;
R = A;    Q = eye(n);
for j = 1:n-2
    w = R(j+1:n, j);
    if norm(w, inf) == 0    break;    end
    omiga = sign(w(1)) * norm(w);    beta = omiga * (omiga + w(1));    w(1) = w(1) + omiga;
    H = eye(n);    H(j+1:n, j+1:n) = H(j+1:n, j+1:n) - w * w' / beta;    R = H * R * H;    Q = Q * H;
end

数值实验

例 3：给定矩阵

用初等反射矩阵正交相似约化后，运行[Q, R] = hessenberg(A)得到Q（正交矩阵）和R（上黑森伯格矩阵）。

>> [Q, R] = hessenberg(A)

4. 上黑森伯格矩阵的 QR 算法

算法原理：设H为不可约上黑森伯格矩阵，μ为的特征值，则 QR 方法中；该算法不能计算复特征值。

算法步骤：

给定为上黑森伯格矩阵，计算

a..

b. 对于k=1,2,⋯,n−1：

(1) 确定旋转变换P(k,k+1)，使

；

(2) 左变换 对于j=k,⋯,n：

.

c. 对于k=1,2,⋯,n−1：

(1) 右变换 对于i=1,2,⋯,k+1：

;

(2).

d..

MATLAB 程序：

function D = D_QR(A, it_max, ep)
% 单步位移QR方法求上黑森伯格矩阵特征值
% it_max为最大迭代次数
[m, n] = size(A);
if m ~= n    error('A非方阵');    end
D = eye(n);
while n > 1
    i = 1;
    while i <= it_max
        if abs(A(n, n-1)) < ep    break;    end
        A = QR_T(A);    i = i + 1;
    end
    D(n, n) = A(n, n);    A = A(1:n-1, 1:n-1);    n = n - 1;
    if n == 1    D(1, 1) = A(1, 1);    end
end
D = diag(D);
function H = QR_T(A)
% 用平面旋转变换做上黑森伯格矩阵的QR分解
n = length(A);    H = A;    st = H(n, n);    H(1, 1) = H(1, 1) - st;
for k = 1:n-1
    H(k+1, k+1) = H(k+1, k+1) - st;
    if abs(H(k+1, k)) <= eps    % 数值稳定的算法
        c(k) = 1; s(k) = 0;
    elseif abs(H(k+1, k)) >= abs(H(k, k))
        t = H(k, k) / H(k+1, k);
        s(k) = 1 / sqrt(1 + t^2); c(k) = s(k) * t;
    else
        t = H(k+1, k) / H(k, k);
        c(k) = 1 / sqrt(1 + t^2); s(k) = c(k) * t;
    end
    H(k:k+1, k:n) = [c(k), s(k); -s(k), c(k)] * H(k:k+1, k:n);
end
for k = 1:n-1
    H(1:k+1, k:k+1) = H(1:k+1, k:k+1) * [c(k), -s(k); s(k), c(k)];
    H(k, k) = H(k, k) + st;
end
H(n, n) = H(n, n) + st;

数值实验

例 4：用 QR 算法求矩阵

的全部特征值：

>> A = [5, -2, -5, -1; 1, 6, 8, 6; 0, 2, 5, -3; 0, 0, 5, 12];    D = D_QR(A, 12, 1e-5)

验证：

>> eig(A)

总结

本文介绍了三种矩阵特征值计算方法：

1. 幂法通过迭代求主特征值及特征向量，可利用原点平移法加速；

2. 反幂法通过求解逆矩阵的特征值来计算最小特征值，同样支持原点平移加速；

3. QR算法将矩阵转化为上黑森伯格形式后迭代求解全部特征值。

文中给出了各算法的MATLAB实现，并通过数值实验验证了有效性，其中幂法迭代次数减少45%，反幂法对接近2的特征值仅需11次迭代，QR算法可精确求得全部特征值。实验表明，特征值分离越好，迭代收敛越快。