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

线性变换 T T T的矩阵 A ∈ P n × n \boldsymbol{A}\in P^{n\times n} A∈Pn×n的对角化,即寻求对角阵 Λ \boldsymbol{\Lambda} Λ,使得 A \boldsymbol{A} A~ Λ \boldsymbol{\Lambda} Λ,需分几步走:
(1)解方程 det ( λ I − A ) = 0 \det(\lambda\boldsymbol{I}-\boldsymbol{A})=0 det(λI−A)=0,得根
λ 1 , λ 2 , ⋯ , λ k ∈ P \lambda_1,\lambda_2,\cdots,\lambda_k\in P λ1,λ2,⋯,λk∈P
为 A \boldsymbol{A} A的特征值;
(2)对每一个特征值 λ i \lambda_i λi,解齐次线性方程组 ( λ i I − A ) x = o (\lambda_i\boldsymbol{I}-\boldsymbol{A})\boldsymbol{x}=\boldsymbol{o} (λiI−A)x=o,得基础解系 α i 1 , α i 2 , ⋯ , α i n i \boldsymbol{\alpha}_{i1},\boldsymbol{\alpha}_{i2},\cdots,\boldsymbol{\alpha}_{in_i} αi1,αi2,⋯,αini, i = 1 , 2 , ⋯ , k i=1,2,\cdots,k i=1,2,⋯,k;
(3)若 n 1 + n 2 + ⋯ + n k = n n_1+n_2+\cdots+n_k=n n1+n2+⋯+nk=n,则
Λ = diag ( λ 1 , ⋯ , λ 1 ⏟ n 1 , ⋯ , λ k , ⋯ , λ k ⏟ n k ) \boldsymbol{\Lambda}=\text{diag}(\underbrace{\lambda_1,\cdots,\lambda_1}_{n_1},\cdots,\underbrace{\lambda_k,\cdots,\lambda_k}_{n_k}) Λ=diag(n1
λ1,⋯,λ1,⋯,nk
λk,⋯,λk)
A \boldsymbol{A} A~ Λ \boldsymbol{\Lambda} Λ。即 T T T在基 α 11 , ⋯ , α 1 n 1 , ⋯ , α k 1 , ⋯ , α k n 1 \boldsymbol{\alpha}_{11},\cdots,\boldsymbol{\alpha}_{1n_1},\cdots,\boldsymbol{\alpha}_{k1},\cdots,\boldsymbol{\alpha}_{kn_1} α11,⋯,α1n1,⋯,αk1,⋯,αkn1下的矩阵为 Λ \boldsymbol{\Lambda} Λ。
numpy.linalg提供了函数eigvals用来计算方阵的特征值,其调用接口为
eigvals(A) \text{eigvals(A)} eigvals(A)
参数A表示方阵 A \boldsymbol{A} A。返回值为 A \boldsymbol{A} A的 n n n个根(包括重根)。需要提起注意的是,此处返回的根有可能是复数。对函数eigvals算出 A \boldsymbol{A} A的每个ℝ中的特征值 λ \lambda λ,调用博文《线性方程组的通解》中定义的mySolve函数,计算 ( λ I − A ) x = o (\lambda\boldsymbol{I}-\boldsymbol{A})\boldsymbol{x}=\boldsymbol{o} (λI−A)x=o的基础解系,记为属于 λ \lambda λ的线性无关的特征向量。把属于各不同特征值的特征向量按序排列,若这些向量的总数恰为 n n n,则构成向量空间的一个基底,且各特征值(包括重复值)构成的对角阵就是 T T T在这个基下的矩阵。
例1 用Python计算矩阵 A = ( 1 2 2 2 1 2 2 2 1 ) ∈ \boldsymbol{A}=\begin{pmatrix}1&2&2\\2&1&2\\2&2&1\end{pmatrix}\in A=
122212221
∈ℝ 3 × 3 ^{3\times3} 3×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行设置矩阵 A \boldsymbol{A} A的数据为A。第8行读取 A \boldsymbol{A} A的阶数为n。第9行调用numpy.linalg的eigvals函数计算 A \boldsymbol{A} A的 n n n个特征值存于w。第10行调用numpy的diag函数用w按升序排列的 n n n个值构造对角阵Lambda。第11行调用numpy的unique函数将w中的 n n n个特征值删掉重复,仅保留不同值。即 det ( λ I − A ) = 0 \det(\lambda\boldsymbol{I}-\boldsymbol{A})=0 det(λI−A)=0的所有重根仅算一个。得到的各不相同的特征值按升序存于v。注意,w中存储的 A \boldsymbol{A} A的特征值都是浮点型的,为能确定重根,需限制精度。np.round(w,decimals=14)调用numpy的round函数将w中的值限制有效位数为14。第12行设置 n n n阶单位阵I,第13行设置 n n n-维零向量o。第14行读取第1个特征值 λ 1 \lambda_1 λ1为lam,第15行调用mySolve函数(见博文《线性方程组的通解》)解齐次线性方程组 ( λ 1 I − A ) = o (\lambda_1\boldsymbol{I}-\boldsymbol{A})=\boldsymbol{o} (λ1I−A)=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]]
即 A = ( 1 2 2 2 1 2 2 2 1 ) ∈ \boldsymbol{A}=\begin{pmatrix}1&2&2\\2&1&2\\2&2&1\end{pmatrix}\in A=
122212221
∈ℝ 3 × 3 ^{3\times3} 3×3在基 α 1 = ( − 1 1 0 ) , α 2 = ( − 1 0 1 ) \boldsymbol{\alpha}_1=\begin{pmatrix}-1\\1\\0\end{pmatrix},\boldsymbol{\alpha}_2=\begin{pmatrix}-1\\0\\1\end{pmatrix} α1=
−110
,α2=
−101
和 α 3 = ( 1 1 1 ) \boldsymbol{\alpha}_3=\begin{pmatrix}1\\1\\1\end{pmatrix} α3=
111
下对角化为 Λ = ( − 1 0 0 0 − 1 0 0 0 5 ) \boldsymbol{\Lambda}=\begin{pmatrix}-1&0&0\\0&-1&0\\0&0&5\end{pmatrix} Λ=
−1000−10005
。
写博不易,敬请支持:
如果阅读本文于您有所获,敬请点赞、评论、收藏,谢谢大家的支持!
代码诚可贵,原理价更高。若为AI学,读正版书好。
更多推荐


所有评论(0)