线性代数Python计算：矩阵对角化

线性变换$T$的矩阵$\mathbf{A}\in P^{n\times n}$的对角化，即寻求对角阵$\mathbf{\Lambda }$，使得$\mathbf{A}$~$\mathbf{\Lambda }$，需分几步走：
(1)解方程$\det (\lambda \mathbf{I}-\mathbf{A})=0$，得根
$${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_k\in P$$
为$\mathbf{A}$的特征值；
(2)对每一个特征值${\lambda }_i$，解齐次线性方程组$({\lambda }_i\mathbf{I}-\mathbf{A})\mathbf{x}=\mathbf{o}$，得基础解系${\mathbf{\alpha }}_{i1},{\mathbf{\alpha }}_{i2},\cdots ,{\mathbf{\alpha }}_{i{n}_i}$，$i=1,2,\cdots ,k$；
(3)若$n_1+n_2+\cdots +n_k=n$，则
$$\mathbf{\Lambda }=\text{diag}(\underset{n_1}{\underbrace{{\lambda }_1,\cdots ,{\lambda }_1}},\cdots ,\underset{n_k}{\underbrace{{\lambda }_k,\cdots ,{\lambda }_k}})$$
$\mathbf{A}$~$\mathbf{\Lambda }$。即$T$在基${\mathbf{\alpha }}_{11},\cdots ,{\mathbf{\alpha }}_{1n_1},\cdots ,{\mathbf{\alpha }}_{k1},\cdots ,{\mathbf{\alpha }}_{k{n}_1}$下的矩阵为$\mathbf{\Lambda }$。
numpy.linalg提供了函数eigvals用来计算方阵的特征值，其调用接口为
$$\text{eigvals(A)}$$
参数A表示方阵$\mathbf{A}$。返回值为$\mathbf{A}$的$n$个根（包括重根）。需要提起注意的是，此处返回的根有可能是复数。对函数eigvals算出$\mathbf{A}$的每个ℝ中的特征值$\lambda$，调用博文《线性方程组的通解》中定义的mySolve函数，计算$(\lambda \mathbf{I}-\mathbf{A})\mathbf{x}=\mathbf{o}$的基础解系，记为属于$\lambda$的线性无关的特征向量。把属于各不同特征值的特征向量按序排列，若这些向量的总数恰为$n$，则构成向量空间的一个基底，且各特征值（包括重复值）构成的对角阵就是$T$在这个基下的矩阵。
例1 用Python计算矩阵$\mathbf{A}=\left(\begin{array}{ccc}1& 2& 2\\ 2& 1& 2\\ 2& 2& 1\end{array}\right)\in$ℝ${}^{3\times 3}$的对角化。

import numpy as np                      			#导入numpy
from fractions import Fraction as F                 #导入Fraction
np.set_printoptions(formatter=                      #设置输出数据格式
                    {'all':lambda x:str(F(x).limit_denominator())})
A=np.array([[1,2,2],                    			#设置矩阵
            [2,1,2],
            [2,2,1]],dtype='float')
n,_=A.shape                             			#A的阶数
w=np.linalg.eigvals(A)                  			#A的特征值
Lambda=np.diag(np.sort(w))              			#对角阵
v=np.unique(np.round(w,decimals=14))    			#特征值唯一化
I=np.eye(n)                             			#单位阵
o=np.zeros((n,1))                       			#零向量
lam=v[0]                                			#第1个特征值
P=(mySolve(lam*I-A,o))[:,1:]            			#属于第1个特征值的特征向量
for lam in v[1:]:                       			#其余个特征值
    X=mySolve(lam*I-A,o)                			#属于特征值的特征向量
    P=np.hstack((P,X[:,1:n]))
print('对角阵：')
print(Lambda)
print('基：')
print(P)

程序的第5~7行设置矩阵$\mathbf{A}$的数据为A。第8行读取$\mathbf{A}$的阶数为n。第9行调用numpy.linalg的eigvals函数计算$\mathbf{A}$的$n$个特征值存于w。第10行调用numpy的diag函数用w按升序排列的$n$个值构造对角阵Lambda。第11行调用numpy的unique函数将w中的$n$个特征值删掉重复，仅保留不同值。即$\det (\lambda \mathbf{I}-\mathbf{A})=0$的所有重根仅算一个。得到的各不相同的特征值按升序存于v。注意，w中存储的$\mathbf{A}$的特征值都是浮点型的，为能确定重根，需限制精度。np.round(w,decimals=14)调用numpy的round函数将w中的值限制有效位数为14。第12行设置$n$阶单位阵I，第13行设置$n$-维零向量o。第14行读取第1个特征值${\lambda }_1$为lam，第15行调用mySolve函数（见博文《线性方程组的通解》）解齐次线性方程组$({\lambda }_1\mathbf{I}-\mathbf{A})=\mathbf{o}$，将所得基础解系（存于返回值的第2列起的各列）置于P中。第16~18行的for循环将其余各特征值所属特征向量依次置于P后，最终得到对角化后的基底。运行程序，输出

对角阵：
[[-1  0 0]
 [ 0 -1 0]
 [ 0  0 5]]
基：
[[-1 -1 1]
 [ 1  0 1]
 [ 0  1 1]]

即$\mathbf{A}=\left(\begin{array}{ccc}1& 2& 2\\ 2& 1& 2\\ 2& 2& 1\end{array}\right)\in$ℝ${}^{3\times 3}$在基${\mathbf{\alpha }}_1=\left(\begin{array}{c}-1\\ 1\\ 0\end{array}\right),{\mathbf{\alpha }}_2=\left(\begin{array}{c}-1\\ 0\\ 1\end{array}\right)$和${\mathbf{\alpha }}_3=\left(\begin{array}{c}1\\ 1\\ 1\end{array}\right)$下对角化为$\mathbf{\Lambda }=\left(\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& 5\end{array}\right)$。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
代码诚可贵，原理价更高。若为AI学，读正版书好。