数值计算方法 Chapter7. 计算矩阵的特征值和特征向量 0. 问题描述1. 幂法 1. 思路2. 规范运算3. 伪代码实现2. 反幂法 1. 思路 & 方法2. 伪代码实现3. 实对称矩阵的Jacobi方法 1. 思路 & 方法2. 伪代码实现0. 问题描述这一章节面对的问题是说，给定一个n 阶矩阵，如何数值求解其特征值，即：
Ax = \lambda x1. 幂法1. 思路幂法的主要思路其实依然还是来源于迭代思想。
显然，对于任意一个向量\vec{x} ，我们总可以将其用n 阶矩阵的一组正交基进行表示，即：
\vec{x} = \sum_{i=1}^{n} x_i \cdot \vec{n_i}其中，\vec{n_i} 为矩阵A 的一个单位向量，有A\vec{n_i} = \lambda_i \vec{n_i} 。令\vec{x}^{(0)} = \sum_{i}x_i \vec{n_i} ，\vec{x}^{(k+1)} = A\vec{x}^{(k)} ，则易有：
\vec{x}^{(k)} = \sum_{i} x_i \cdot \lambda_i^k \cdot \vec{n_i}对应的模长为：
||\vec{x}^{(k)}|| = \sqrt{\sum_i x_i^2 \lambda_i^{2k}}我们考虑当k \to \infty 时\vec{x}^{(k)} 的方向，不妨设max(|\lambda_i|) = \lambda 。
不妨设\lambda = |\lambda_1| \geq |\lambda_2| \geq ... \geq |\lambda_n| ，则我们可以分情况讨论有：
假设有且仅有一个正的\lambda_i 使得|\lambda_i| = \lambda \lim_{k \to \infty} \frac{\vec{x}^{(k)}}{||\vec{x}^{(k)}||} = \vec{n_1}
此时有：
\vec{x}^{(k+1)} = A\vec{x}^{(k)} = \lambda \vec{x}^{(k)} 假设满足|\lambda_i| = \lambda 条件的i 共有m 个，且均有\lambda_i > 0 \lim_{k \to \infty} \frac{\vec{x}^{(k)}}{||\vec{x}^{(k)}||} = \sum_{i}^{m}\frac{x_i}{\sqrt{\sum_i^m x_i^2}}\vec{n_i}
此时同样有：
\vec{x}^{(k+1)} = A\vec{x}^{(k)} = \lambda \vec{x}^{(k)} 假设有m 个正的|\lambda_i| = \lambda ，l 个负的|\lambda_i| = \lambda \lim_{k \to \infty} \frac{\vec{x}^{(k)}}{||\vec{x}^{(k)}||} = \sum_{i}^{m}\frac{x_i}{\sqrt{\sum_i^{m+l} x_i^2}}\vec{n_i} + \sum_{j}^{l}\frac{(-1)^{k} \cdot x_j}{\sqrt{\sum_i^{m+l} x_j^2}}\vec{n_j}
此时，虽然\vec{x}^{(k+1)} / \vec{x}^{(k)} 不稳定，但是依然有：
\vec{x}^{(k+2)} = AA\vec{x}^{(k)} = \lambda^2 \vec{x}^{(k)}不过需要注意的是，这里讨论的只是一般的情况，其基于的假设是说所有的x_i \neq 0 ，如果恰好存在某些分量上的投影为0， 即某些x_i = 0 ，那么上述讨论会发生一定的变化甚至失效。
而且，如上述分析，通过幂法，我们只能够获得一般矩阵当中绝对值最大的一个特征值\lambda ，无法获取其所有的特征值，这个也需要注意一下。
2. 规范运算基于上述思路，我们给出幂法计算的规范运算方法：
\left\{ \begin{aligned} \vec{y}^{(k)} &= \vec{x}^{(k)} / ||\vec{x}^{(k)}|| \\ \vec{x}^{(k+1)} &= A \cdot \vec{y}^{(k)} \end{aligned} \right.通过上述迭代，我们即可求取矩阵A 的绝对值最大的特征值。
3. 伪代码实现给出python伪代码实现如下：
代码语言：javascript复制def power_fn(A, epsilon=1e-6):
n = len(A)
x = [1 for _ in range(n)]
for _ in range(10**3):
s = math.sqrt(sum(t*t for t in x))
y = [t/s for t in x]
z = [0 for _ in range(n)]
for i in range(n):
for j in range(n):
z[i] += A[i][j] * y[j]
lamda = sum([a/b for a, b in zip(z, y)]) / n
delta = max(abs(a-b) for a, b in zip(z, x))
x = z
if delta < epsilon:
break
return lamda当然，这里只考虑了仅有单一符号的最大绝对值\lambda 的情况，如果存在绝对值相同且符号相反的最大特征值\lambda ，上述代码会失效，需要进行一些复杂操作还有判断，这里就不过多展开了。
2. 反幂法1. 思路 & 方法反幂法的思路和幂法其实大差不差，不过幂法是直接正向的进行迭代，即：
\vec{x}^{(k+1)} = A \vec{x}^{(k)}反幂法的迭代方式则是：
\vec{x}^{(k+1)} = A^{-1} \vec{x}^{(k)}或者说：
\vec{x}^{(k)} = A \vec{x}^{(k+1)}解法而言，我们只需要在上述幂法的基础上叠加上一章节当中多元方程组求解的方法即可。
需要额外说明的是，由于这里使用的迭代与之前的幂法是相反的，因此，这里求解的是A^{-1}
当中绝对值最大的特征值，也就是A 当中绝对值最小的特征值。
而同样的，这里的额外隐性条件就是需要矩阵A 是满秩的，否则矩阵不存在逆矩阵，上述方程\vec{x}^{(k)} = A \vec{x}^{(k+1)} 可能无解。
2. 伪代码实现同样的，这里我们给出最简单情况下（即A 满秩且且仅存在一个最小的绝对值特征根的情况），反幂法的python伪代码实现：
代码语言：javascript复制def rev_power_fn(A, epsilon=1e-6):
n = len(A)
x = [1 for _ in range(n)]
for _ in range(10**3):
s = math.sqrt(sum(t*t for t in x))
y = [t/s for t in x]
z = loose_iter(A, y)
lamda = sum([a/b for a, b in zip(y, z)]) / n
delta = max(abs(a-b) for a, b in zip(z, x))
x = z
if delta < epsilon:
break
return lamda3. 实对称矩阵的Jacobi方法1. 思路 & 方法如前所述，幂法和反幂法本质上都是通过迭代的思路找一个稳定的特征向量，然后通过特征向量来求特征值。
因此，他们只能求取矩阵的某一个特征值，无法对矩阵的全部特征值进行求解。如果要对矩阵的全部特征值进行求解，上述方法就会失效。
但是，对于一些特殊的矩阵，即实对称矩阵，事实上我们是可以对其全部的特征值进行求解的，一种典型的方法就是Jacobi方法。
本质上来说，Jacobi方法依然还是进行迭代，不过其迭代的思路则是不断地对矩阵进行酉变换，使之收敛到一个对角矩阵上面，此时对角矩阵的各个对角元就是原矩阵的特征值。
具体而言，已知Ax = \lambda x ，有正交矩阵U 满足UU^{T} = I ，则有：
UAU^{T} \cdot Ux = \lambda \cdot Ux此时，\lambda 依然是矩阵UAU^{T} 的特征值。
如果经过一系列变换之后：
U_kU_{k-1}...U_{1}AU_1^{T}U_2^{T}...U_{k}^T = diag(\lambda_1, ..., \lambda_n)则\lambda_1, ..., \lambda_n 即为矩阵A 的全部特征值。
剩下的问题就是如何求解这些矩阵U ，Jacobi方法给出的一种可行的思路是通过Givens矩阵，即：
G(p, q, \theta) = \begin{pmatrix} 1 \\ & ... \\ & & cos\theta & ... & sin\theta \\ & & ... & & ... \\ & & -sin\theta & ... & cos\theta \\ & & & & & ... \\ & & & & & & 1 \end{pmatrix}也就是说，只对p,q 两行的元素进行角度为\theta 的旋转变换。
令:
B = G(p,q,\theta) \cdot A \cdot G^{T}(p,q,\theta)则有：
\left\{ \begin{aligned} b_{i, j} &= a_{i, j} \\ b_{i,p} &= b_{p,i} = a_{p, i}cos\theta - a_{q, i}sin\theta \\ b_{i,q} &= b_{q,i} = a_{p, i}sin\theta + a_{q, i}cos\theta \\ b_{p,p} &= a_{p,p}cos^{2}\theta + a_{q,q}sin^{2}\theta - a_{p,q}sin2\theta \\ b_{q,q} &= a_{p,p}sin^{2}\theta + a_{q,q}cos^{2}\theta + a_{p,q}sin2\theta \\ b_{p,q} &= b_{q,p} = a_{p,q} cos2\theta + \frac{a_{p,p} - a_{q,q}}{2}sin2\theta \end{aligned} \right.我们取合适的\theta ，使得b_{q,p} = 0 ，则有：
\left\{ \begin{aligned} \sum_{i\neq j} b_{i,j}^2 &= \sum_{i\neq j} a_{i,j}^2 - 2a_{p,q}^2 \\ \sum_{i} b_{i,i}^2 &= \sum_{i} a_{i,i}^2 + 2a_{p,q}^2 \end{aligned} \right.可以看到，非对角元的元素绝对值会越来越小。
因此，经过足够次数的迭代，可以将原始矩阵A 变换成为一个特征值相同的近对角矩阵。
而为了进一步提升迭代的速度，可以优先选择绝对值最大的非对角元进行迭代消去。
2. 伪代码实现同样的，我们给出python伪代码实现如下：
代码语言：javascript复制def jacobi_eigen_fn(A, epsilon=1e-6):
n = len(A)
assert(len(A[0]) == n)
assert(A[i][j] == A[j][i] for i in range(n) for j in range(i+1, n))
finished = False
for _ in range(1000):
for p in range(n):
for q in range(p+1, n):
if A[p][p] != A[q][q]:
theta = 0.5 * math.atan(2*A[p][q] / (A[q][q] - A[p][p]))
else:
theta = math.pi / 4
ap = deepcopy(A[p])
aq = deepcopy(A[q])
for i in range(n):
if i == p or i == q:
continue
A[i][p] = A[p][i] = ap[i] * math.cos(theta) - aq[i] * math.sin(theta)
A[i][q] = A[q][i] = ap[i] * math.sin(theta) + aq[i] * math.cos(theta)
A[p][p] = ap[p] * math.cos(theta)**2 + aq[q] * math.sin(theta)**2 - ap[q] * math.sin(2*theta)
A[q][q] = ap[p] * math.sin(theta)**2 + aq[q] * math.cos(theta)**2 + ap[q] * math.sin(2*theta)
A[p][q] = A[q][p] = ap[q] * math.cos(2 * theta) + (ap[p] - aq[q])/2 * math.sin(2*theta)
finished = all(A[i][j] < epsilon for i in range(n) for j in range(i+1, n))
if finished:
break
if finished:
break
if finished:
break
return [A[i][i] for i in range(n)]