数值计算方法 Chapter7. 计算矩阵的特征值和特征向量

2022-06-20 09:00:38 浏览数 (1)
  • 数值计算方法 Chapter7. 计算矩阵的特征值和特征向量
    • 0. 问题描述
    • 1. 幂法
      • 1. 思路
      • 2. 规范运算
      • 3. 伪代码实现
    • 2. 反幂法
      • 1. 思路 & 方法
      • 2. 伪代码实现
    • 3. 实对称矩阵的Jacobi方法
      • 1. 思路 & 方法
      • 2. 伪代码实现

0. 问题描述

这一章节面对的问题是说，给定一个n 阶矩阵，如何数值求解其特征值，即：

Ax = lambda x

1. 幂法

1. 思路

幂法的主要思路其实依然还是来源于迭代思想。

显然，对于任意一个向量vec{x} ，我们总可以将其用n 阶矩阵的一组正交基进行表示，即：

vec{x} = sum_{i=1}^{n} x_i cdot vec{n_i}

其中，vec{n_i} 为矩阵A 的一个单位向量，有Avec{n_i} = lambda_i vec{n_i} 。令vec{x}^{(0)} = sum_{i}x_i vec{n_i} ，vec{x}^{(k 1)} = Avec{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| ，则我们可以分情况讨论有：

  1. 假设有且仅有一个正的lambda_i 使得|lambda_i| = lambda
lim_{k to infty} frac{vec{x}^{(k)}}{||vec{x}^{(k)}||} = vec{n_1}

此时有：

vec{x}^{(k 1)} = Avec{x}^{(k)} = lambda vec{x}^{(k)}
  1. 假设满足|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)} = Avec{x}^{(k)} = lambda vec{x}^{(k)}
  1. 假设有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)} = AAvec{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 lamda

3. 实对称矩阵的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 \ & ... \ & & costheta & ... & sintheta \ & & ... & & ... \ & & -sintheta & ... & costheta \ & & & & & ... \ & & & & & & 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}costheta - a_{q, i}sintheta \ b_{i,q} &= b_{q,i} = a_{p, i}sintheta a_{q, i}costheta \ b_{p,p} &= a_{p,p}cos^{2}theta a_{q,q}sin^{2}theta - a_{p,q}sin2theta \ b_{q,q} &= a_{p,p}sin^{2}theta a_{q,q}cos^{2}theta a_{p,q}sin2theta \ b_{p,q} &= b_{q,p} = a_{p,q} cos2theta frac{a_{p,p} - a_{q,q}}{2}sin2theta end{aligned} right.

我们取合适的theta ，使得b_{q,p} = 0 ，则有：

left{ begin{aligned} sum_{ineq j} b_{i,j}^2 &= sum_{ineq 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)]
对象存储 serverless

0 人点赞