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

codename_cys
239 阅读6分钟

本文已参与「新人创作礼」活动，一起开启掘金创作之路。

0. 问题描述

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

Ax=λxAx = \lambda x

1. 幂法

1. 思路

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

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

x⃗=∑i=1nxi⋅ni⃗\vec{x} = \sum_{i=1}^{n} x_i \cdot \vec{n_i}

其中，ni⃗\vec{n_i}为矩阵AA的一个单位向量，有Ani⃗=λini⃗A\vec{n_i} = \lambda_i \vec{n_i}。

令x⃗(0)=∑ixini⃗\vec{x}^{(0)} = \sum_{i}x_i \vec{n_i}，x⃗(k+1)=Ax⃗(k)\vec{x}^{(k+1)} = A\vec{x}^{(k)}，则易有：

x⃗(k)=∑ixi⋅λik⋅ni⃗\vec{x}^{(k)} = \sum_{i} x_i \cdot \lambda_i^k \cdot \vec{n_i}

对应的模长为：

∣∣x⃗(k)∣∣=∑ixi2λi2k||\vec{x}^{(k)}|| = \sqrt{\sum_i x_i^2 \lambda_i^{2k}}

我们考虑当k→∞k \to \infty时x⃗(k)\vec{x}^{(k)}的方向，不妨设max(∣λi∣)=λmax(|\lambda_i|) = \lambda。

不妨设λ=∣λ1∣≥∣λ2∣≥...≥∣λn∣\lambda = |\lambda_1| \geq |\lambda_2| \geq ... \geq |\lambda_n|，则我们可以分情况讨论有：

  1. 假设有且仅有一个正的λi\lambda_i使得∣λi∣=λ|\lambda_i| = \lambda

    lim⁡k→∞x⃗(k)∣∣x⃗(k)∣∣=n1⃗\lim_{k \to \infty} \frac{\vec{x}^{(k)}}{||\vec{x}^{(k)}||} = \vec{n_1}

    此时有：

    x⃗(k+1)=Ax⃗(k)=λx⃗(k)\vec{x}^{(k+1)} = A\vec{x}^{(k)} = \lambda \vec{x}^{(k)}

  2. 假设满足∣λi∣=λ|\lambda_i| = \lambda条件的ii共有mm个，且均有λi>0\lambda_i > 0

    lim⁡k→∞x⃗(k)∣∣x⃗(k)∣∣=∑imxi∑imxi2ni⃗\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}

    此时同样有：

    x⃗(k+1)=Ax⃗(k)=λx⃗(k)\vec{x}^{(k+1)} = A\vec{x}^{(k)} = \lambda \vec{x}^{(k)}

  3. 假设有mm个正的∣λi∣=λ|\lambda_i| = \lambda，ll个负的∣λi∣=λ|\lambda_i| = \lambda

    lim⁡k→∞x⃗(k)∣∣x⃗(k)∣∣=∑imxi∑im+lxi2ni⃗+∑jl(−1)k⋅xj∑im+lxj2nj⃗\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}

    此时，虽然x⃗(k+1)/x⃗(k)\vec{x}^{(k+1)} / \vec{x}^{(k)}不稳定，但是依然有：

    x⃗(k+2)=AAx⃗(k)=λ2x⃗(k)\vec{x}^{(k+2)} = AA\vec{x}^{(k)} = \lambda^2 \vec{x}^{(k)}

不过需要注意的是，这里讨论的只是一般的情况，其基于的假设是说所有的xi≠0x_i \neq 0，如果恰好存在某些分量上的投影为0， 即某些xi=0x_i = 0，那么上述讨论会发生一定的变化甚至失效。

而且，如上述分析，通过幂法，我们只能够获得一般矩阵当中绝对值最大的一个特征值λ\lambda，无法获取其所有的特征值，这个也需要注意一下。

2. 规范运算

基于上述思路，我们给出幂法计算的规范运算方法：

{y⃗(k)=x⃗(k)/∣∣x⃗(k)∣∣x⃗(k+1)=A⋅y⃗(k)\left\{ \begin{aligned} \vec{y}^{(k)} &= \vec{x}^{(k)} / ||\vec{x}^{(k)}|| \\ \vec{x}^{(k+1)} &= A \cdot \vec{y}^{(k)} \end{aligned} \right.

通过上述迭代，我们即可求取矩阵AA的绝对值最大的特征值。

3. 伪代码实现

给出python伪代码实现如下：

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. 思路 & 方法

反幂法的思路和幂法其实大差不差，不过幂法是直接正向的进行迭代，即：

x⃗(k+1)=Ax⃗(k)\vec{x}^{(k+1)} = A \vec{x}^{(k)}

反幂法的迭代方式则是：

x⃗(k+1)=A−1x⃗(k)\vec{x}^{(k+1)} = A^{-1} \vec{x}^{(k)}

或者说：

x⃗(k)=Ax⃗(k+1)\vec{x}^{(k)} = A \vec{x}^{(k+1)}

解法而言，我们只需要在上述幂法的基础上叠加上一章节当中多元方程组求解的方法即可。

需要额外说明的是，由于这里使用的迭代与之前的幂法是相反的，因此，这里求解的是A−1A^{-1}当中绝对值最大的特征值，也就是AA当中绝对值最小的特征值。

而同样的，这里的额外隐性条件就是需要矩阵AA是满秩的，否则矩阵不存在逆矩阵，上述方程x⃗(k)=Ax⃗(k+1)\vec{x}^{(k)} = A \vec{x}^{(k+1)}可能无解。

2. 伪代码实现

同样的，这里我们给出最简单情况下（即AA满秩且且仅存在一个最小的绝对值特征根的情况），反幂法的python伪代码实现：

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=λxAx = \lambda x，有正交矩阵UU满足UUT=IUU^{T} = I，则有：

UAUT⋅Ux=λ⋅UxUAU^{T} \cdot Ux = \lambda \cdot Ux

此时，λ\lambda依然是矩阵UAUTUAU^{T}的特征值。

如果经过一系列变换之后：

UkUk−1...U1AU1TU2T...UkT=diag(λ1,...,λn)U_kU_{k-1}...U_{1}AU_1^{T}U_2^{T}...U_{k}^T = diag(\lambda_1, ..., \lambda_n)

则λ1,...,λn\lambda_1, ..., \lambda_n即为矩阵AA的全部特征值。

剩下的问题就是如何求解这些矩阵UU，Jacobi方法给出的一种可行的思路是通过Givens矩阵，即：

G(p,q,θ)=(1...cosθ...sinθ......−sinθ...cosθ...1)G(p, q, \theta) = \begin{pmatrix} 1 \\ & ... \\ & & cos\theta & ... & sin\theta \\ & & ... & & ... \\ & & -sin\theta & ... & cos\theta \\ & & & & & ... \\ & & & & & & 1 \end{pmatrix}

也就是说，只对p,qp,q两行的元素进行角度为θ\theta的旋转变换。

令:

B=G(p,q,θ)⋅A⋅GT(p,q,θ)B = G(p,q,\theta) \cdot A \cdot G^{T}(p,q,\theta)

则有：

{bi,j=ai,jbi,p=bp,i=ap,icosθ−aq,isinθbi,q=bq,i=ap,isinθ+aq,icosθbp,p=ap,pcos2θ+aq,qsin2θ−ap,qsin2θbq,q=ap,psin2θ+aq,qcos2θ+ap,qsin2θbp,q=bq,p=ap,qcos2θ+ap,p−aq,q2sin2θ\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，使得bq,p=0b_{q,p} = 0，则有：

{∑i≠jbi,j2=∑i≠jai,j2−2ap,q2∑ibi,i2=∑iai,i2+2ap,q2\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.

可以看到，非对角元的元素绝对值会越来越小。

因此，经过足够次数的迭代，可以将原始矩阵AA变换成为一个特征值相同的近对角矩阵。

而为了进一步提升迭代的速度，可以优先选择绝对值最大的非对角元进行迭代消去。

2. 伪代码实现

同样的，我们给出python伪代码实现如下：

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)]
avatar
文章
阅读
粉丝