第二十四篇 特征值问题编程基础：特征值特征向量的求解和性质

特征值方程

在分析结构稳定性或振动系统的固有频率时，经常会出现这种情况。我们必须找到一个向量{x}，当它与[a]相乘时，得到它自身的标量倍。这个λ倍数称为[A]的“特征值”，我们将看到n阶矩阵有n个这样的λ。从物理上讲，它们可能代表振荡的频率。还有与每个特征值λ相关的n个向量{x}。这些被称为“特征向量”。从物理上讲，它们可能代表振荡的模态振型。下面是一个具体的例子：

重新排列后可以写成下面的形式

这组线性联立方程只有在系数的行列式为零时才可能有非0解

行列式展开得

这叫做“特征多项式”。显然，解特征值方程的一种方法是将其简化为一个n次特征多项式，并使用前一章的方法来求其根。有时这样做，也许作为整个解过程的一部分，但它本身通常不是求解特征值方程的最佳方法。
上面得特征多项式有一个简单得因式分解

所以特征值解为4，1，-8
对于任意矩阵[A]，特征多项式都可能产生虚根和实根。所以要把问题限制在具有实数特征值的矩阵上，物理约束通常意味着矩阵是对称的和“正定的”，在这种情况下，所有的特征值都是实数和正数。找到一个特征值后，其相关的特征向量可以通过求解一组线性联立方程来找到。例如，将λ = 1代入方程得到

经过高斯消元得第一个阶段

由上面可以看出，行列式得值为0，方程数比未知数数量少，可见方程组是线性相关的。然而，我们可以从上面的方程第二行和第三行看出，x2与x3的比值是2:1，通过代换方程的第一行，x1与x2的比值是1:1。所以当特征值λ = 1时，任意一个特征向量的比值x1: x2: x3 = 2:2:1。类似的操作也可以用来求两个特征值对应的特征向量。

特征向量的正交化和规范化

特征向量有一个独特的方向，但它们的大小是任意的。两边同时乘以一个任意常数α，得到。

在该修正方程中，α{x}和{x}是数组[A]中对应特征值λ的等效特征向量。两种常用的特征向量正交化方法。对于下面的对称矩阵

具有特征值

对应的特征向量

特征向量正交化的一种简便形式包括缩放长度，或者称为“L2范数”或“欧几里德范数”。因此对于一个特征向量的n个分量，我们像下面这个缩放这个向量

它是累加全部分量的平方和再求跟得到的
因此，对于第一个特征向量，将除以下面这个长度值

对于其他两个特征向量，得到正交特征向量为

对称矩阵的特征向量互相具有“正交性”。也就是说，来自同一个矩阵的任意两个不同特征向量的点积等于零，而特征向量与自身的点积为其长度的平方。因此，将欧几里德范数设为单位1的一个便利的特点是是特征向量点积要么等于零，要么等于单位1。

特征向量正交化的一种简单的替代形式是除以具有最大(绝对)幅度的分量。得到了一个最大分量等于单位1的正交化特征向量。对于之前给出的特征向量，第一个特征向量除以√2，第二个特征向量不变，因第三个特征向量将被-√2整除得到正交化向量。

特征值和特征向量的性质

本部分所描述的数值方法利用了很多矩阵变换对特征值和特征向量影响的定理。下面总结了这些关系中使用较多的一些。

如果一个n × n矩阵[A]具有实数特征值λ1， λ2，···，λn， [I]为单位矩阵，p为标量位移，则

n × n矩阵对角线的和称为“迹”，也等于特征值的和，因此

一个n×n矩阵的特征值的乘积等于这个矩阵的行列式，因此

一个n × n三角形矩阵的对角线是特征值。

给定两个方阵[A]和[B]。

[A][B]和[B][A]将有相同的特征值。
如果{x}是[A][B]的特征向量，则[B]{x}是[B][A]的特征向量。
解释：

然后两边都乘以B得到

给定两个方阵[A]和[P]，其中[P]是非0的。

[A]和[P]−1[A][P]具有相同的特征值。
如果{x}是[A]的特征向量，那么[P]−1{x}是[P]−1[A][P]的特征向量。
解释：

通过两边都乘以[P]−1，得到

嵌入矩阵形式[P][P]-1=[I]，得到

如果[A]是一个n×n矩阵，其特征值和特征向量为λ1{x1}， λ2{x2}，···，λn{xn}。[P]是一个矩阵，其列是[a]的使其欧几里德范数等于单位的特征向量缩放，那么

一 QR原理

理论依据：任意一个非奇异矩阵（满秩的方阵）A都可以分解为一个正交矩阵Q和一个上三角矩阵R的乘积，且当R对角元符号确定时，分解是唯一的。QR分解是一种迭代方法，迭代格式如下：
当Ak基本收敛到为上三角矩阵时，迭代完成，此时主对角元素就是特征值。

特别地：当A是对称阵的时候，Ak是对角阵Λ，Q=Qk-1Qk-2…Q1就是其正交特征向量矩,有Q^TAQ=Ak=Λ，即A正交对角化与Ak。

如何理解？我们看下图公式:

所以，QR迭代过程从数学的角度来想其实就是不断正交化的过程。

二 QR算法步骤

1.Householder变换进行QR分解

反射矩阵：任取单位向量w,反射矩阵H=E-2WW^T ,显然HH^T =E，H是正交阵

定理：任取两个模长相等的的向量x,y，一定存在一个反射矩阵H，使得Hx=y，
此时w=(x-y)/(|x-y|)（向量的差除以向量差的模）

应用：现在我们取矩阵的一列为x,m=|x|,y=m*[1,0,0,…0]^T 根据上面的定理求出H，使得Hx=y,是不是通过正交变化就把那一列化成了[m,0,0,0]^T ,这样就达到了将下三角元素全化为0的效果。看下图，举个例子来说明QR分解过程：

看懂上述过程就知道，Householder变换是利用了反射定理，经过n-1轮正交变换，将下三角元素全部化为0，从而得到上三角矩阵R，将所有H矩阵左乘运算再转置得到正交矩阵Q，即A=QR

我们看看QR分解的代码：

#QR分解
def qrSplit(A):
    n=A.shape[0]#A的维度
    Q=[[]]
    R=A
    for i in range(0,n-1):
        B=R
        if i!=0:
            #删除一行一列,得n-1阶子阵
            B=B[i:,i:]
        #取第一列向量
        x=B[:,0]
        #向量摸长
        m=np.linalg.norm(x)
        #生成一个模长为m，其余项为0的向量y
        y=[0 for j in range(0,n-i)]
        y[0]=m
        #计算householder反射矩阵
        #w = (x-y)/||x-y||
        w=x-y
        w=w/np.linalg.norm(w)
        #H=E-2*WT*W
        H=np.eye(n-i)-2*np.dot(w.reshape(n-i,1),w.reshape(1,n-i))#H是个正交矩阵
        #第一次计算不需对正交正H升维
        if i==0: 
            #第一次迭代
            Q=H
            R=np.dot(H,R)
        else:
            #因为降了维度，所以要拼接单位矩阵
            D=np.c_[np.eye(i),np.zeros((i,n-i))]
            H=np.c_[np.zeros((n-i,i)),H]
            H=np.r_[D,H]
            #迭代计算正交矩阵Q和上三角R
            Q=np.dot(H,Q)
            R=np.dot(H,R)
    Q=Q.T
    return [Q,R]

我们测试一下QR分解后是否能还原：

A=np.array([1,2,3,4,2,1,2,3,3,2,1,2,4,3,2,1])
A=A.reshape(4,4)
print('A原来的样子')
print(A)
qr = qrSplit(A)
print('打印Q,R')
print(qr[0])
print(qr[1])
print('打印Q*R')
print(np.dot(qr[0],qr[1]))

很明显，R下三角元素都非常小，可以认为是0.

2.QR迭代

上一步完成了QR分解，QR迭代就非常简单了，看代码：

#QR迭代求特征值特征向量
def  qrEgis(A):
    # QR迭代(尽量让它多迭代几次，以至于AK收敛为上三角)
    qr = []
    n = A.shape[0]  # A的维度
    Q = np.eye(n)
    for i in range(0, 100):
        # A=QR
        qr = qrSplit(A)
        # 将Q右边边累成
        Q = np.dot(Q,qr[0])
        # A1=RQ
        A = np.dot(qr[1], qr[0])

    AK = np.dot(qr[0], qr[1])
    #把e取出来
    e=[]
    for i in range(0,n):
        e.append(AK[i][i])
    #对特征值按升序排序，特征向量与其对应
    for i in range(0,n-1):
        min=e[i]
        for j in range(i+1,n):
            if e[j]<min:
                min=e[j]
                #交换特征值
                tmp=e[i]
                e[i]=e[j]
                e[j]=tmp
                #交换特征向量
                r=np.copy(Q[:,i])
                Q[:,i]=Q[:,j]
                Q[:,j]=r;
    return [e,Q]

我们输入一个对称阵，同时去求出它的特征值与特征向量，测试一下。
并且与numpy自带的求解特征值、特征向量的做个对比，看代码：

A=np.array([1,2,3,4,2,1,2,3,3,2,1,2,4,3,2,1])
A=A.reshape(4,4)
egis =qrEgis(A)
print('自己写的QR分解')
print(egis[0]);
print('......')
print(egis[1]);
print('numpy自带的分解器')
e,u=np.linalg.eigh(A);
print(e)
print('.....')
print(u)

可以看出自己写的QR分解法和numpy自带的分解器求出的特征值是一样的。
特征向量在数值上一样，符号不一样并没有关系，因为X和-X都是A的特征向量。

注意：如果A不是对称阵，QR法只能求出全部特征值，不能同时求出特征向量。
但是，已知特征值求特征向量可以采用反幂法。
关于反幂法，学过数值分析的知道，其实也是一个迭代过程。

引言

• 参考链接1
• 参考链接2
• 参考链接3
• 参考链接4

Jacobi法求特征值特征向量

求解使用对象：实对称矩阵$R$,性质有$R{R}^{-1}=I$ 即是 $R^T=R^{-1}$。
　　求解思路：实对称矩阵性质与旋转矩阵$R$性质相同，故看做是将旋转矩阵转化为三个自由度。变换过程也是利用旋转矩阵相乘来达到。化为对角线性的过程就是旋转的过程，对角线上的元素就是特征值，用来变换的旋转矩阵的列向量就是特征向量。
　　用处蛮多的：对于协方差矩阵、海塞矩阵的近似$J^{T}J$、旋转矩阵的分解均可以使用此方法，精度高，速度快，过程不涉及求逆。分解后的形式：$A=V\Lambda V^T$，没有求逆，在后续计算也比较方便。
　　名字为什么叫$Jacobi$法没有深入。