一、幂法和反幂法
1、幂法
由于求解特征方程$|\boldsymbol{I}-\lambda\boldsymbol{A}|=0$的计算量巨大，所以在求解特征值时，要从特征值的定义
$$ \boldsymbol{Ax}=\lambda\boldsymbol{x} $$
入手。设 $n \times n$ 阶实矩阵 $A$ 的特征值 $\lambda_i \ (i=1,2,\cdots,n)$ 满足：
$$ |\lambda_1| > |\lambda_2| \geq \cdots \geq |\lambda_n| $$
且与 $\lambda_i\ (i=1,2,\cdots,n)$ 相应的特征向量 $\boldsymbol{u}_1, \boldsymbol{u}_2, \cdots, \boldsymbol{u}_n$ 线性无关。给定初始向量 $\boldsymbol{x}^{(0)} \neq 0$，由迭代公式：
$$ \boldsymbol{x}^{(k+1)} = \boldsymbol{A} \boldsymbol{x}^{(k)} \qquad (k=0,1,2,\cdots) $$
产生向量序列 $\{\boldsymbol{x}^{(k)}\}$，当 $k$ 充分大时，有 $\lambda_1 \approx x^{(k+1)}_i / x^{(k)}_i$，相应的特征向量为 $\boldsymbol{x}^{(k+1)}$。
证明当$k$充分大时，上述结论成立：
设 $\|\boldsymbol{u}_i\| = 1 \ (i=1,2,\cdots,n)$，且 $\boldsymbol{u}_i$ 线性无关，故必存在 $n$ 个不全为零的数 $\alpha_i\ (i=1,2,\cdots,n)$，使得：
$$ \boldsymbol{x}^{(0)} = \sum_{i=1}^{n} \alpha_i \boldsymbol{u}_i $$
得：
$$ \begin{aligned} \boldsymbol{x}^{(k+1)} &= A^{k+1} \boldsymbol{x}^{(0)} = A^{k+1} \left(\sum_{i=1}^n \alpha_i \boldsymbol{u}_i \right) = \sum_{i=1}^n \alpha_i A^{k+1} \boldsymbol{u}_i \\ &= \sum_{i=1}^n \alpha_i \lambda_i^{k+1} \boldsymbol{u}_i \end{aligned} $$
即：
$$ \boldsymbol{x}^{(k+1)} = \lambda_1^{k+1} \left[\alpha_1 \boldsymbol{u}_1 + \left(\frac{\lambda_2}{\lambda_1}\right)^{k+1} \alpha_2 \boldsymbol{u}_2 + \cdots + \left(\frac{\lambda_n}{\lambda_1}\right)^{k+1} \alpha_n \boldsymbol{u}_n \right] $$
设 $\alpha_1 \neq 0$, 由 $|\lambda_1| > |\lambda_i|, \ (i = 2, 3, \cdots, n)$ 得：
$$ \lim_{k \to \infty} \left( \frac{\lambda_i}{\lambda_1} \right)^{k+1} \alpha_i \boldsymbol{u}_i = 0 $$
于是：
$$ \lim_{k \to \infty} \sum_{i=2}^n \left( \frac{\lambda_i}{\lambda_1} \right)^{k+1} \alpha_i \boldsymbol{u}_i = 0 $$
故只要 $k$ 充分大, 就有：
$$ \boldsymbol{x}^{(k+1)} = \lambda_1^{k+1} \left[ \alpha_1 \boldsymbol{u}_1 + \sum_{i=2}^n \left( \frac{\lambda_i}{\lambda_1} \right)^{k+1} \alpha_i \boldsymbol{u}_i \right] \approx \lambda_1^{k+1} \alpha_1 \boldsymbol{u}_1 $$
因此, 可把 $\boldsymbol{x}^{(k+1)}$ 作为与 $\lambda_1$ 相应的特征向量的近似。由公式：
$$ \boldsymbol{x}^{(k+1)} \approx \lambda_1^{k+1} \alpha_1 \boldsymbol{u}_1, \qquad \boldsymbol{x}^{(k)} \approx \lambda_1^{k} \alpha_1 \boldsymbol{u}_1 $$
即可得到：
$$ \lambda_1 \approx \frac{x^{(k+1)}_i}{x^{(k)}_i} \qquad (i = 1, 2, \cdots, n) $$
利用上述结论，可以不断迭代，得到：
$$ \boldsymbol{x}^{k} = \boldsymbol{A}^k \boldsymbol{x}^{(0)} \qquad (k=0,1,2,\cdots) $$
注意到其中$\boldsymbol{A}^k$在$k$较大的情况下会出现两种状况：
- 如果$\boldsymbol{A}$的范数大于$1$，则会导致计算过程中溢出。
- 如果$\boldsymbol{A}$的范数小于$1$，则会导致在计算过程中变成$0$。
根据上面的推导过程，以及矩阵2-范数的定义，为了避免这种情形，实际计算时每次迭代所求得的向量都要归一化，因此实际计算使用公式为：
$$ \begin{cases} \boldsymbol{y}^{(k)} = \dfrac{\boldsymbol{x}^{(k)}}{x^{(k)}_r}, \qquad \left| x^{(k)}_r \right| = \max\limits_{1 \leq i \leq n} \left| x^{(k)}_i \right| \qquad (k = 0,1,2,\cdots) \\ \boldsymbol{x}^{(k+1)} = A \boldsymbol{y}^{(k)} \\ \lambda_1 \approx x^{(k+1)}_r \end{cases} $$
当我们计算$|\boldsymbol{Ax}^{(k+1)}-\lambda\boldsymbol{x}^{(k+1)}|$的结果小于判定值时，结束循环，判定为其已经收敛，得到特征向量$\boldsymbol{x}^{(k+1)}$，在此之后，一般不做除法来获得特征值，而是使用瑞利商：
$$ \lambda_1=\frac{\boldsymbol{x}^\mathrm{T}\boldsymbol{Ax}}{\boldsymbol{x}^\mathrm{T}\boldsymbol{x}} $$
这样所得的结果更加稳定。
2、反幂法
在实际应用是，我们需要寻找一个最接近于数$m$的特征值，可以通过求解$\boldsymbol{A}-m\boldsymbol{E}$的最小特征值（绝对值最小），从而将所得结果加$m$后获得答案。
设 $\boldsymbol{A}$ 为 $n \times n$ 阶非奇异矩阵，$\lambda, \boldsymbol{u}$ 为 $\boldsymbol{A}$ 的特征值与相应的特征向量，则 $\boldsymbol{A}^{-1}$ 的特征值是 $\boldsymbol{A}$ 的特征值的倒数，而相应的特征向量不变，即
$$ \boldsymbol{A}^{-1}\boldsymbol{u} = \frac{1}{\lambda}\boldsymbol{u} $$
因此，若对矩阵 $\boldsymbol{A}^{-1}$ 用幂法，即可计算出 $\boldsymbol{A}^{-1}$ 的按模最大的特征值，其倒数恰为 $\boldsymbol{A}$ 的按模最小的特征值。这就是反幂法的基本思想。
由于$\boldsymbol{A}^{-1}$的计算比较麻烦，而且往往不能保持矩阵$\boldsymbol{A}$的一些好的性质（如稀疏性），因此反幂法在实际运算时以求解方程组：
$$ \boldsymbol{Ax}^{(k+1)}=\boldsymbol{x}^{(k)} $$
代替，同时在每次迭代的过程中做归一化处理。
3、代码
下面是使用反幂法求解特定特征值（距离给出 lambda0 最近的特征值）的MATLAB代码：
function [lambda, v] = solveEigenpair(A, lambda0)
n = size(A, 1);
B = A - lambda0 * eye(n);
[L, U, P] = lu(B);
v = rand(n, 1);
v = v / norm(v);
tol = 1e-12;
while true
v_new = U \ (L \ (P * v));
v_new = v_new / norm(v_new);
lambda = v_new' * A * v_new;
if norm(A * v_new - lambda * v_new) < tol
break;
end
v = v_new;
end
end二、雅可比（Jacobi）方法
Jacobi 方法是求实对称矩阵的全部特征值及相应特征向量的一种方法 。它基于以下两个数学结论：
- 正交相似变换：任意实对称矩阵 $\boldsymbol{A}$ 均存在正交矩阵 $\boldsymbol{Q}$，使得 $\boldsymbol{Q}^T\boldsymbol{AQ} = \mathrm{diag}(\lambda_1, \lambda_2, \cdots, \lambda_n)$，其中 $\lambda_i$ 为特征值，$\boldsymbol{Q}$ 的列向量为特征向量。
- Frobenius 范数不变性：在正交变换下，矩阵元素的平方和保持不变 。设 $\boldsymbol{B} = \boldsymbol{Q}^T\boldsymbol{AQ}$，则 $\sum a_{ij}^2 = \sum b_{ij}^2$。
该方法的核心思想是通过一系列正交变换（旋转变换），不断减小非对角元素的平方和，使矩阵趋向于对角阵。下面介绍旋转变换矩阵 $V_{ij}(\varphi)$，旋转变换矩阵是将单位阵的第 $i, j$ 行/列交点处的元素替换为三角函数值得到：
$$ \boldsymbol{V}_{ij}(\varphi) = \begin{bmatrix} 1 & & & & & \\ & \ddots & & & & \\ & & \cos \varphi & \cdots & \sin \varphi & \\ & & \vdots & 1 & \vdots & \\ & & -\sin \varphi & \cdots & \cos \varphi & \\ & & & & & \ddots & \\ & & & & & & 1 \end{bmatrix} \begin{matrix} \\ \\ i \\ \\ j \\ \\ \\ \end{matrix} $$
选取旋转角 $\varphi$ 使得变换后 $a_{ij}^{(1)} = 0$。若 $a_{ii} \neq a_{jj}$，则取：
$$ \tan 2\varphi = \frac{2a_{ij}}{a_{ii} - a_{jj}} \quad (|\varphi| \leq \frac{\pi}{4}) $$
雅各比方法运算量大，且不能保持矩阵的特殊形状，所以该方法只适合求解中小型稠密矩阵。
三、QR 方法
1、基本 QR 算法
QR 算法是计算中小型矩阵全部特征值最有效的方法之一 。其基本思想是利用矩阵的 QR 分解进行相似变换迭代，可以证明，任意非奇异实矩阵都可以分解成一个正交矩阵$\boldsymbol{Q}$和一个上三角矩阵$\boldsymbol{R}$的乘积，并且当$\boldsymbol{R}$的对角元符号确定时，分解时唯一的，其分解格式为：
$$ \begin{cases} \boldsymbol{A}^{(k)} = \boldsymbol{Q}_k \boldsymbol{R}_k \\ \boldsymbol{A}^{(k+1)} = \boldsymbol{R}_k \boldsymbol{Q}_k \end{cases} \quad (k=1, 2, \cdots) $$
其中 $\boldsymbol{Q}_k$ 为正交阵，$\boldsymbol{R}_k$ 为上三角阵 。迭代产生的序列 $\{\boldsymbol{A}^{(k)}\}$ 相似于原矩阵，并最终收敛于上三角形（或分块上三角形）矩阵。因为上三角阵的主对角元是该矩阵的特征值，故$k$充分大时，$\boldsymbol{A}^{(k)}$的主对角元可以看作是原矩阵的特征值近似。
使用施密特（Schmidt）正交化方法对 $n$ 阶非奇异实矩阵 $\boldsymbol{A}$ 进行 QR 分解的通式推导过程如下：
设 $\boldsymbol{A}$ 为 $n$ 阶非奇异实矩阵，记为列向量形式：
$$ \boldsymbol{A} = [\boldsymbol{a}_1, \boldsymbol{a}_2, \cdots, \boldsymbol{a}_n] $$
其中 $\boldsymbol{a}_j = (a_{1j}, a_{2j}, \cdots, a_{nj})^T \quad (j=1, 2, \cdots, n)$。利用施密特正交化构造正交向量组 $\boldsymbol{b}_1, \boldsymbol{b}_2, \cdots, \boldsymbol{b}_n$ ，对于第一列向量 $\boldsymbol{a}_1$，直接进行单位化：
$$ \boldsymbol{b}_1 = \boldsymbol{a}_1 / \|\boldsymbol{a}_1\| $$
对于后续列向量 $\boldsymbol{a}_k(k=2, 3, \cdots, n)$，先求正交向量 $\boldsymbol{b}_k'$，从 $\boldsymbol{a}_k$ 中减去其在已求出的正交向量 $\boldsymbol{b}_1, \cdots, \boldsymbol{b}_{k-1}$上的投影：
$$ \boldsymbol{b}_k' = \boldsymbol{a}_k - \sum_{i=1}^{k-1} \langle \boldsymbol{a}_k, \boldsymbol{b}_i \rangle \boldsymbol{b}_i $$
得到长度为 1 的正交向量 $\boldsymbol{b}_k$：
$$ \boldsymbol{b}_k = \boldsymbol{b}_k' / \|\boldsymbol{b}_k'\| $$
将上述迭代公式改写为 $\boldsymbol{a}_k$的线性组合形式：
$$ \boldsymbol{a}_k = \langle \boldsymbol{a}_k, \boldsymbol{b}_1 \rangle \boldsymbol{b}_1 + \cdots + \langle \boldsymbol{a}_k, \boldsymbol{b}_{k-1} \rangle \boldsymbol{b}_{k-1} + \|\boldsymbol{b}_k'\| [cite_start]\boldsymbol{b}_k $$
根据该线性组合关系，矩阵 $\boldsymbol{A}$ 可以写成正交矩阵 $\boldsymbol{Q}$ 与上三角矩阵 $\boldsymbol{R}$ 的乘积：
$$ \boldsymbol{A} = [\boldsymbol{b}_1, \boldsymbol{b}_2, \cdots, \boldsymbol{b}_n] \begin{bmatrix} \|\boldsymbol{a}_1\| & \langle \boldsymbol{a}_2, \boldsymbol{b}_1 \rangle & \cdots & \langle \boldsymbol{a}_n, \boldsymbol{b}_1 \rangle \\ & \|\boldsymbol{b}_2'\| & \cdots & \langle \boldsymbol{a}_n, \boldsymbol{b}_2 \rangle \\ & & \ddots & \vdots \\ & & & \|\boldsymbol{b}_n'\| \end{bmatrix} = \boldsymbol{QR} $$
2、Householder 变换
为了减少运算量，通常先通过 Householder 变换将一般矩阵化为拟上三角形（Hessenberg 矩阵）。设 $\|\boldsymbol{w}\|_2 = 1$，则称 $\boldsymbol{H} = \boldsymbol{I} - 2\boldsymbol{ww}^T$ 为 Householder 矩阵。$\boldsymbol{H}$ 是实对称的正交矩阵（$\boldsymbol{H} = \boldsymbol{H}^T = \boldsymbol{H}^{-1}$），且 $\det(\boldsymbol{H}) = -1$。它可以将任一非零向量 $\boldsymbol{x}$ 变换为单位向量的倍数（即实现坐标轴旋转/反射），用于矩阵消零。
该矩阵存在一个重要性质：设 $\boldsymbol{x}, \boldsymbol{y}$ 为 $\mathbb{R}^n$ 中任意两个非零向量，且 $\|\boldsymbol{y}\|_2 = 1$ ，则存在 Householder 矩阵 $\boldsymbol{H}$，使得：
$$ \boldsymbol{H}\boldsymbol{x} = \pm \|\boldsymbol{x}\|_2 \boldsymbol{y} $$
该定理表明，对于任一非零向量 $\boldsymbol{x}$，都可以构造一个 Householder 变换，将其变成事先给定的单位向量的倍数。在数值计算中，通常取 $\boldsymbol{y} = \boldsymbol{e}_i$（单位坐标向量），使 $\boldsymbol{x}$ 变换后只有一个分量不为零。
为了避免在 $\boldsymbol{x}$ 与 $\boldsymbol{e}_i$ 方向接近时产生较大的舍入误差，实际计算中总取以下通式构造 $\boldsymbol{w}$：
$$ \boldsymbol{w} = \frac{\boldsymbol{x} + \|\boldsymbol{x}\|_2 \boldsymbol{e}_i \text{sign}(x_i)}{\|\boldsymbol{x} + \|\boldsymbol{x}\|_2 \boldsymbol{e}_i \text{sign}(x_i)\|_2} $$
3、化一般矩阵为拟上三角矩阵
拟上三角形矩阵（上海森堡阵）指主对角线下方的元素从第二行开始全为零的矩阵 。其形式如下：
$$ \boldsymbol{H} = \begin{bmatrix} h_{11} & h_{12} & \cdots & h_{1n} \\ h_{21} & h_{22} & \cdots & h_{2n} \\ 0 & h_{32} & \cdots & h_{3n} \\ \vdots & \ddots & \ddots & \vdots \\ 0 & \cdots & h_{n,n-1} & h_{nn} \end{bmatrix} $$
首先构造第一个 Householder 矩阵 $\boldsymbol{H}_1$，选取 $\boldsymbol{A}$ 矩阵第一列对角线以下的元素组成向量 $\boldsymbol{a}_1 = (a_{21}, a_{31}, \cdots, a_{n1})^T$ 。根据定理，构造一个 $n-1$ 阶的 Householder 矩阵 $\tilde{\boldsymbol{H}}_1$，使得变换后的向量 $\tilde{\boldsymbol{H}}_1 a_1$ 只有第一个分量不为零，即 $\tilde{\boldsymbol{H}}_1 \boldsymbol{a}_1 = \alpha(1, 0, \cdots, 0)^T$ 。将 $\tilde{\boldsymbol{H}}_1$ 嵌入到 $n$ 阶单位阵中：
$$ \boldsymbol{H}_1 = \begin{bmatrix} 1 & 0 \\ 0 & \tilde{H}_1 \end{bmatrix} $$
计算 $\boldsymbol{H}_1 \boldsymbol{A} \boldsymbol{H}_1$，此时所得矩阵的第一列除前两个元素外，其余元素均为零。对于 $k = 2, 3, \cdots, n-2$，依次构造 $\boldsymbol{H}_k$。$\boldsymbol{H}_k$ 的形式为：
$$ \boldsymbol{H}_k = \begin{bmatrix} I_k & 0 \\ 0 & \tilde{H}_k \end{bmatrix} $$
其中 $\boldsymbol{I}_k$ 是 $k$ 阶单位阵，$\tilde{\boldsymbol{H}}_k$ 是对当前矩阵第 $k$ 列对角线以下元素进行处理的 Householder 矩阵。经过 $n-2$ 次正交相似变换后，得到拟上三角形矩阵 $H$：
$$ \boldsymbol{H} = \boldsymbol{H}_{n-2} \cdots \boldsymbol{H}_2 \boldsymbol{H}_1 \boldsymbol{A} \boldsymbol{H}_1 \boldsymbol{H}_2 \cdots \boldsymbol{H}_{n-2} $$
4、拟上三角矩阵的QR分解
由于拟上三角形矩阵 $\boldsymbol{H}$ 的特殊形状（仅在次对角线上有非零元），通常使用 $n-1$ 个旋转变换（即 Givens 变换）将其化为上三角形矩阵，从而实现 QR 分解 。通过一系列旋转矩阵 $\boldsymbol{V}_{i, i-1}$ 依次消去拟上三角矩阵 $\boldsymbol{H}$ 的次对角线元素 $h_{21}, h_{32}, \cdots, h_{n, n-1}$。
取旋转矩阵 $\boldsymbol{V}_{21}$ 消去 $h_{21}$ ：
$$ \boldsymbol{V}_{21}=\begin{bmatrix}\cos\varphi_1&\sin\varphi_1&0&\cdots&0\\-\sin\varphi_1&\cos\varphi_1&0&\cdots&0\\&&1&\cdots&0\\&&&\ddots&\vdots\\&&&&1\end{bmatrix} $$
计算参数：$r_1 = \sqrt{h_{11}^2 + h_{21}^2}$，$cos \varphi_1 = h_{11}/r_1$，$sin \varphi_1 = h_{21}/r_1$。变换$H^{(2)} = V_{21}H$，此时 $h_{21}^{(2)} = 0$ 。
循环操作，假设前 $k-1$ 步已完成，取 $V_{k+1, k}$ 消去 $H^{(k)}$ 的次对角线元 $h_{k+1, k}^{(k)}$ 。计算参数$r_k = \sqrt{(h_{kk}^{(k)})^2 + (h_{k+1, k}^{(k)})^2}$，$\cos \varphi_k = h_{kk}^{(k)}/r_k$，$\sin \varphi_k = h_{k+1, k}^{(k)}/r_k$ 。变换$\boldsymbol{H}^{(k+1)} = \boldsymbol{V}_{k+1, k}\boldsymbol{H}^{(k)}$ 。经过最多 $n-1$ 次旋转变换后，矩阵变为上三角矩阵 $R$ ：
$$ \boldsymbol{R} = \boldsymbol{H}^{(n)} = \boldsymbol{V}_{n, n-1}\boldsymbol{V}_{n-1, n-2} \cdots \boldsymbol{V}_{21}\boldsymbol{H} $$
由于旋转矩阵均为正交矩阵，根据正交阵的性质：
$$ \boldsymbol{H} = \boldsymbol{V}_{21}^T \boldsymbol{V}_{32}^T \cdots \boldsymbol{V}_{n, n-1}^T \boldsymbol{R} = \boldsymbol{QR} $$
其中，$\boldsymbol{Q} = \boldsymbol{V}_{21}^T \boldsymbol{V}_{32}^T \cdots \boldsymbol{V}_{n, n-1}^T$ 仍为正交矩阵。完成这一过程的运算量约为$4n^2$，比一般矩阵的QR分交少了一个数量级。
5、代码
下面为使用 Householder 矩阵优化后的QR分解的MATLAB代码：
function [Q, H] = hessenberg(A)
n = size(A, 1);
Q = eye(n);
H = A;
for i = 1:n-2
x = H(i+1:n, i);
e = zeros(n-i, 1);
e(1) = 1;
s = sign(x(1));
if s == 0
s = 1;
end
w = (x + norm(x, 2) * e * s) / norm(x + norm(x, 2) * e * s, 2);
H0 = eye(n-i) - 2 * (w * w');
H1 = eye(n);
H1(i+1:n, i+1:n) = H0;
Q = Q * H1;
H = H1 * H * H1;
end
end
function [Q, R] = getQR(H)
n = size(H, 1);
Q = eye(n);
R = H;
for i = 1:n-1
a = R(i, i);
b = R(i+1, i);
if b == 0
c = 1;
s = 0;
else
r = hypot(a, b);
c = a / r;
s = b / r;
end
ri = R(i, i:n);
rj = R(i+1, i:n);
R(i, i:n) = c * ri + s * rj;
R(i+1, i:n) = -s * ri + c * rj;
qi = Q(:, i);
qj = Q(:, i+1);
Q(:, i) = c * qi + s * qj;
Q(:, i+1) = -s * qi + c * qj;
end
end
function values = eigen(A)
n = size(A, 1);
[~, H] = hessenberg(A);
tol = 1e-5;
while norm(H - diag(diag(H)), 'fro') >= tol
mu = H(n, n);
B = H(1:n, 1:n) - mu * eye(n);
[Q, R] = getQR(B);
H = R * Q + mu * eye(n);
end
values = real(diag(H));
end